Signaling Events Downstream of AHR Activation That Contribute to Toxic Responses: The Functional Role of an AHR-Dependent Long Noncoding RNA (slincR) Using the Zebrafish Model

Background: A structurally diverse group of chemicals, including dioxins [e.g., 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD)] and polycyclic aromatic hydrocarbons (PAHs), can xenobiotically activate the aryl hydrocarbon receptor (AHR) and contribute to adverse health effects in humans and wildlife. In the zebrafish model, repression of sox9b has a causal role in several AHR-mediated toxic responses, including craniofacial cartilage malformations; however, the mechanism of sox9b repression remains unknown. We previously identified a long noncoding RNA, sox9b long intergenic noncoding RNA (slincR), which is increased (in an AHR-dependent manner) by multiple AHR ligands and is required for the AHR-activated repression of sox9b. Objective: Using the zebrafish model, we aimed to enhance our understanding of the signaling events downstream of AHR activation that contribute to toxic responses by identifying: a) whether slincR is enriched on the sox9b locus, b) slincR’s functional contributions to TCDD-induced toxicity, c) PAHs that increase slincR expression, and d) mammalian orthologs of slincR. Methods: We used capture hybridization analysis of RNA targets (CHART), qRT-PCR, RNA sequencing, morphometric analysis of cartilage structures, and hemorrhaging screens. Results: The slincR transcript was enriched at the 5′ untranslated region (UTR) of the sox9b locus. Transcriptome profiling and human ortholog analyses identified processes related to skeletal and cartilage development unique to TCDD-exposed controls, and angiogenesis and vasculature development unique to TCDD-exposed zebrafish that were injected with a splice-blocking morpholino targeting slincR. In comparison to TCDD exposed control morphants, slincR morphants exposed to TCDD resulted in abnormal cartilage structures and a smaller percentage of animals displaying the hemorrhaging phenotype. In addition, slincR expression was significantly increased in six out of the sixteen PAHs we screened. Conclusion: Our study establishes that in zebrafish, slincR is recruited to the sox9b 5′ UTR to repress transcription, can regulate cartilage development, has a causal role in the TCDD-induced hemorrhaging phenotype, and is up-regulated by multiple environmentally relevant PAHs. These findings have important implications for understanding the ligand-specific mechanisms of AHR-mediated toxicity. https://doi.org/10.1289/EHP3281


Introduction
Recent advances in genetics and molecular biotechnology have resulted in a shift in the field of toxicology from observations of apical end points to investigations of the molecular initiating events in the adverse outcome pathway (National Research Council 2007). The aryl hydrocarbon receptor (AHR) is a conserved ligand-dependent transcription factor that is a member of the basic helix-loop-helix-PER-ARNT-SIM protein family. The AHR is known to be xenobiotically activated by a structurally diverse group of chemicals, including dioxins and polycyclic aromatic hydrocarbons (PAHs; Denison and Nagy 2003;Hahn et al. 2017). Upon ligand binding, the AHR translocates to the nucleus and heterodimerizes with the aryl hydrocarbon nuclear translocator (ARNT), inducing ligand-specific transcriptional changes (Beischlag et al. 2008). In vertebrates, the AHR also plays key roles in regulating a variety of physiological responses, such as immune and developmental processes, suggesting that too little or too much AHR activity can lead to adverse health effects (Esser and Rannug 2015;Murray et al. 2014). In zebrafish, the activation of the AHR by 2,3,7,8-tetrachlorodibenzo-p-dioxin (TCDD) and numerous PAHs results in the dysregulation of hundreds of genes and is associated with lethality, wasting, hemorrhaging, and defects in the skeletal and cardiovascular systems (Chlebowski et al. 2017;Denison and Nagy 2003;Garcia et al. 2018;Goodale et al. 2015).
The best-characterized high-affinity ligand used to study AHRdependent toxicity is TCDD, a halogenated aromatic hydrocarbon. TCDD and related halogenated aromatic hydrocarbons are extremely toxic to most vertebrate species. The developmental zebrafish is among the most sensitive vertebrate models, making it ideal for investigating AHR signal transduction and function (Tanguay et al. 1999). Ahr-knockout mice and zebrafish mutants have demonstrated that the AHR is required for TCDD-induced toxicity (Fernandez-Salguero et al. 1995;Garcia et al. 2018;Goodale et al. 2012). In zebrafish, Ahr2 and Arnt1 are the functional orthologs of mammalian AHR and ARNT (Antkiewicz et al. 2006;Carney et al. 2006). Mammals and zebrafish share many toxic responses due to their well-conserved genomes, cell types, tissues, and organ systems Howe et al. 2013). Although the majority of downstream genes responsible for specific phenotypes of AHR-dependent toxicity remain largely unknown, the zebrafish model has been used to identify repression of sox9b as a key element responsible for producing craniofacial malformations, and to a lesser extent, cardiovascular malformations (Hofsteen et al. 2013;Xiong et al. 2008). Cardiovascular and craniofacial malformations are hallmarks of TCDD-induced developmental toxicity in fish and mammals (Couture et al. 1990;Henry et al. 1997). SOX9 (human ortholog of zebrafish Sox9b) is a conserved transcription factor that acts as the master regulator of cartilage development (Akiyama et al. 2002).
The mechanism by which AHR activation leads to the repression of sox9b and subsequent developmental toxicities is still unknown; however, we have previously reported the identification of a long noncoding RNA (sox9b long intergenic noncoding RNA, slincR) that is up-regulated by multiple Ahr2 ligands and is located adjacent to sox9b . Long noncoding RNAs (lncRNAs) are defined as transcripts equal to or greater than 200 base pairs (bps) long that do not appear to code for proteins. LncRNAs regulate an array of biological processes, including development, response to stress, embryonic viability, and cancer progression (Gutschner and Diederichs 2012;Li and Chang 2014;Nakagawa 2016). We have used two independent ahr2-mutant zebrafish lines (ahr2 hu3335 and ahr2 osu1 ) to demonstrate that ahr2 is required for the TCDD-induced increase in slincR expression Garcia et al. 2018). Furthermore, ahr2-null zebrafish exposed to TCDD do not display a significant increase in slincR; slincR is required for typical tissue-specific expression of sox9b during normal development and is expressed in tissues with sox9b-essential functions, such as the jaw and snout regions, eyes, and brain . Importantly, slincR expression is also required for the TCDD-induced repression of sox9b.
PAHs can activate the AHR and are ubiquitous environmental contaminants that are primarily created during the incomplete combustion of organic material (e.g., wood, coal, petrol, and oil; Shen et al. 2013). PAHs are of major concern due to their prevalence and potential toxic effects on ecosystems and human health (Perera 1997). PAHs can be genotoxic and/or carcinogenic, and in utero exposure has been associated with adverse birth outcomes (Boström et al. 2002;Dejmek et al. 2000;Perera et al. 2005). In human cohorts from urban environments, prenatal exposure to PAHs has also been associated with increased asthma prevalence and impaired cognitive development in offspring (Edwards et al. 2010;Tang et al. 2012).
The present study is broken down into the following four parts: "PART I: Transcriptional Regulation by slincR," "PART II: slincR Contributions to TCDD Toxicity," "PART III: Investigation of PAHs That Increase slincR Expression," and "PART IV: Identification of Potential slincR Mouse and Human Orthologs." In "PART I: Transcriptional Regulation by slincR," we investigated the mechanism of slincR repression of sox9b using capture hybridization analysis of RNA targets (CHART) to determine if slincR is enriched at the sox9b locus. We hypothesized that slincR binds (indirectly or directly) to the sox9b locus to repress sox9b expression, which contributes to the TCDD-induced cartilage malformation phenotype. We used a previously published sox9b morpholino to determine whether sox9b expression regulates the expression of slincR. We also analyzed the concentration-response relationship between developmental TCDD exposure and the gene expression of sox9b and slincR. In "PART II: slincR Contributions to TCDD Toxicity," to investigate slincR's role in the TCDD-induced toxicity pathway, we performed RNA sequencing (RNA-seq) and Gene Ontology (GO) term enrichment analyses on 48-h post fertilization (hpf) control and slincR morphants exposed to 0.1% DMSO or 1 ng=mL TCDD. We hypothesized that with the reduction/absence of slincR expression (slincR morphants), there will be a relief in sox9b repression, and cartilage-associated GO processes will fail to be significantly enriched. To investigate the phenotypic impact of slincR expression on TCDD-induced jaw malformations, we measured the cartilage of 72-hpf control and slincR morphants treated with 0.1% DMSO or 1 ng=mL TCDD as described in Xiong et al. (2008). Due to the GO term enrichment of angiogenesis and vasculature-related processes in the slincR morphants, we performed a hemorrhaging screen on TCDD-exposed zebrafish at 48 hpf. For both the cartilage and hemorrhaging phenotypic analyses, we hypothesized that the slincR morphants exposed to TCDD will have a less severe toxicity phenotype and/or will be significantly different in comparison with the control morphant datasets. In "PART III: Investigation of PAHs That Increase slincR Expression," we determined if slincR expression was affected in response to sixteen different environmentally relevant PAHs. Finally, in "PART IV: Identification of Potential slincR Mouse and Human Orthologs," we mined multiple RNA-seq datasets to identify the potential mouse and human orthologs of slincR.

Materials and Methods
All gene-specific primers and probes were obtained from Integrated DNA Technologies and are listed in Supplemental Material, Table S1.

Fish Husbandry
The Tropical 5D line of zebrafish (Danio rerio) were purchased from 5D Tropical, Inc. The animals were reared according to Institutional Animal Care and Use Committee protocols at the Sinnhuber Aquatic Research Laboratory, Oregon State University. Adult fish were raised in a recirculating water system (28 ± 1 C) with a 14-h:10-h light:dark schedule. Adult fish were fed GEMMA Micro 300 or 500 (Skretting, Inc.) twice a day. Larval and juvenile fish were fed GEMMA Micro 75 and 150, respectively, thrice a day (Barton et al. 2016). Spawning and embryo collection were conducted as described in Westerfield (2007). Briefly, zebrafish were housed in densities of ∼ 500 fish per 50-gallon tank in recirculating water supplemented with Instant Ocean salts at 28°C. The embryos were collected using a spawning funnel, staged, and maintained in an incubator at 28°C.

Waterborne Exposure
Shield-stage ( ∼ 6 hpf) embryos were exposed to 0, 0.0625, 0.125, 0.25, 0.5, or 1 ng=mL TCDD (311 nM, 95.3% purity; SUPELCO Solutions Within) or vehicle (0.1% DMSO) with gentle rocking for 1 h in 20-mL glass vials (10 embryos=mL). The 1-h duration of TCDD exposure was previously shown sufficient to produce the full TCDD-induced toxicity phenotypes (Henry et al. 1997). The 1 ng=mL concentration was selected to match our previously published TCDD exposure paradigm and results in 99-100% of 120-hpf zebrafish displaying the expected TCDD-induced toxicity phenotype (e.g., cartilage and heart malformations and edema) and an approximate log 2 ðfold changeÞ of −1 in sox9b expression Mathew et al. 2008). We know that TCDD partitions quickly into the embryonic compartment, and the half-life of TCDD is quite long due to a lack of metabolism; thus, this 1 ng=mL TCDD exposure paradigm is sufficient to fully activate the AHR (DeVito and Birnbaum 1995). During the exposures, vials were also gently inverted every 15 min to ensure proper mixing. Embryos were rinsed three times with embryo media and then raised in 100-mm Petri dishes (at 28°C) in approximately 50 mL of embryo media until collection at 48 hpf or 72 hpf, depending on the experimental assay. Embryo media consisted of 15 mM NaCl, 0:5 mM KCl, 1 mM MgSO 4 , 0:15 mM KH 2 PO 4 , 0:05 mM Na 2 HPO 4 , and 0:7 mM NaHCO 3 (Westerfield 2000). The number of embryos raised in a Petri dish ranged from 10-50, depending on the subsequent assay.
PART I: Transcriptional Regulation by slincR qRT-PCR CHART. CHART is a method to identify lncRNA binding sites in the genome (Simon 2013). Zebrafish embryo preparation for CHART was adapted from Bogdanović et al. (2013), and CHART protocol and buffers are from Simon (2013). A detailed protocol is provided (Supplemental Material, S1: CHART Sample Preparation). In brief, zebrafish embryos were exposed to 0.1% DMSO or 1 ng=mL TCDD, as described above. At 24 hpf, embryos were dechorionated with pronase (31:77 lg=lL). At 48 hpf, embryos were screened for malformations, and the morphologically normal animals were euthanized by overdose of buffered tricaine methanesulfonate  by prolonged immersion and monitored under a dissecting microscope until their hearts stopped beating (approximately 15 min). Whole 48-hpf zebrafish embryos were fixed with 4% PFA for 15 min at room temperature, and nuclei were isolated using a Dounce homogenizer and centrifugation. The isolated nuclei were crosslinked using 3.2% PFA and sonicated for 10 min in a Bioruptor ® Pico with cycles of 30 s on and 30 s off. Further, 108 pmol of a 60-bp biotin-labeled capture probe targeting exon 1 of slincR was added to 100 lL of the sonication solution and incubated overnight at room temperature on an end-overend rotator. The control probe comprised of the sense sequence of the same region, and the input samples did not contain a hybridization probe. The capture probes were washed four times with wash buffer and captured using Dynabeads™ MyOne™ Streptavidin T1 (Thermo Fisher Scientific, cat. no. 65601). The input sample was set aside during the washing phase. The samples were de-cross-linked using proteinase K (Ambion™, Thermo Fisher Scientific, cat. no. AM2548) and incubated at 55°C for 1 h, followed by 1 h and 40 min at 65°C. The enriched DNA and RNA were isolated using the ZYMO ZR-Duet™ DNA/RNA miniprep kit (cat. no. D7001) according to manufacturer's instructions. To determine slincR enrichment, 2 lL of RNA pull down was converted to cDNA using a SuperScript™ VILO™ cDNA synthesis kit (Thermo Fisher Scientific, cat. no. 11754050) according to manufacturer's instructions. The genomic DNA pull down was used to determine if slincR is enriched at the sox9b locus. Both the cDNA and genomic DNA enrichment were diluted at 1:20 and were analyzed with primers targeting the slincR transcript and sox9b locus via qRT-PCR using SYBR ® green PCR master mix (Thermo Fisher Scientific, cat. no. 4312704) according to manufacturer's instructions.
Data analysis. Each condition had 3 biological replicates, where 1 replicate consisted of approximately 500 48-hpf zebrafish. The qRT-PCR data were first normalized to 1% input control, such that 6.644 cycles (i.e. dilution factor log 2 ð100Þ) was subtracted from the cycle threshold value (CT; i.e., number of PCR replication cycles required for the sample signal to exceed background levels) of the diluted input and used to calculate the DCT for the two probe sets (DCT = CT½probe − CT½1% input-6:644). For DNA fold enrichment, we next adjusted relative to the sense probe (DDCT = DCT½slincR-probe − DCT½sense-probe), and then fold enrichment was calculated (2 −DDCT ). The RNA yield was calculated using the following equation (2 −DDCT × 100%). We assigned samples that did not amplify (no enrichment) a CT value of 40. The data were graphed using GraphPad Prism 7.02 software. The products were also analyzed via 1.2% agarose gel electrophoresis.
RNA extraction and mRNA quantification. Total RNA was extracted from 48-hpf whole embryos using RNAzol ® (Molecular Research Center, Inc.) and a bullet blender with 0:5 mM zirconium oxide beads, (Next Advance) as recommended by the manufacturer. The RNA was purified using the Direct-zol MiniPrep kit (Zymo Research) and included an in-column DNase 1 digestion. RNA quality and quantity were assessed using a BioTek ® Synergy™ Mx microplate reader with the Gen5™ Take3™ module.
For the sox9b-MO experiment, total RNA (500 ng) was reverse transcribed into cDNA with random primers using the ABI High-Capacity cDNA Reverse Transcription Kit (Thermo Fisher). In addition, qRT-PCR was performed using a StepOnePlus™ Real-Time PCR System (Applied Biosystems). The 20 lL reactions consisted of 10 lL 2X SYBR ® Green Master Mix (Applied Biosystems), 0:4 lL each of 10 lM forward and reverse primers, and 15 ng cDNA. Expression values were normalized to b-actin and analyzed with the 2 −DDCT method as described in Livak and Schmittgen (2001). We used the following calculation to determine morpholino knockdown efficiencies: Statistical analysis. Each biological sample consisted of RNA from 20 pooled 48-hpf zebrafish with 4 biological replicates per condition (n = 4). Results were statistically analyzed and graphed with GraphPad Prism (version 7), and significance was determined using a one-way ANOVA with Tukey post hoc test or Kruskal-Wallis test with Dunn's post hoc test for data that passed or failed normality using the Shapiro-Wilk test, respectively.
For the concentration-response experiment, zebrafish embryos were exposed to six concentrations of TCDD (0, 0.0625, 0.125, 0.25, 0.5, and 1:0 ng=mL) at the shield stage as described in the "Waterborne Exposure" section of the "Methods." The 10-lL one-step qRT-PCR reactions were set up consisting of 5 lL SYBR ® Green Master Mix and 0:08 lL reverse transcriptase enzyme mix (Power SYBR ® Green RNA-to-CT™ 1-Step Kit; Applied Biosystems), 0:2 lL each of 10 lM forward and reverse primers, and 15 ng RNA per reaction. The QuantStudio 5 Real-Time PCR System (Thermo Fisher Scientific) was used under the following cycling conditions: reverse transcription at 48°C for 30 min, denaturation and activation of SYBR ® polymerase at 95°C for 10 min, followed by 40 cycles of amplification (95°C for 15 s, 60°C for 1 min). Expression values were normalized to b-actin and analyzed with the 2 −DDCT method as described in Livak and Schmittgen (2001). The 0 ng=mL TCDD concentration served as the calibrator.
Statistical analysis. Each sample consisted of RNA from 20 pooled 48-hpf zebrafish with 4 biological replicates per condition. Results were statistically analyzed using R (version 3.4.1) in the RStudio (version 1.0.143) integrated development environment (IDE; R version 3.4.1; R Development Core Team; RStudio Team 2016). The data were log 2 -transformed and initially assessed for equal variance and normality using Levene's test and the Shapiro-Wilk test, respectively. Statistical significance was determined using a one-way ANOVA with a Dunnett post hoc test (n = 4, p < 0:001) from the multcomp package (version 1.4-8), and the data were graphed using ggplot2 (version 3.0.0; Hothorn et al. 2008;Wickham 2016).
Evaluation of mortality and 17 physical endpoints at 120 hpf . Zebrafish embryos were exposed to six concentrations of TCDD (0, 0.0625, 0.125, 0.25, 0.5, 1:0 ng=mL) at the shield stage as described in the "Waterborne Exposure" section of the "Methods." The animals were loaded into round-bottom 96-well plates, with one embryo in 100 lL embryo media per well. We measured mortality and a suite of 17 end points (yolk sac edema, body axis, eye, snout, jaw, otic vesicle, pericardial edema, brain, somite, pectoral fin, caudal fin, circulation, pigmentation, trunk length, swim bladder, notochord distortion, and alterations in touch response) at 120 hpf as described in Truong et al. (2014).
Statistical analysis. For each concentration, n = 32 across two 96-well plates (16 animals per plate per treatment group). The 0 ng=mL TCDD concentration served as the control. Due to the low category counts, the data were analyzed using a Fisher's exact test as it does not make distributional assumptions. To control for the family-wise error rate, we applied the Bonferroni correction for multiple comparisons (p < 0:01).
120-hpf larval photomotor response (LPR) assay. Zebrafish embryos were exposed to 0, 0.0625, or 0:125 ng=mL of TCDD at the shield stage ( ∼ 6 hpf), as described in the "Waterborne exposure" section of the Methods. At 120 hpf, only embryos that were phenotypically normal underwent a larval photomotor response assay. The larval photomotor response was conducted using the ViewPoint Zebrabox system (ViewPoint Behavior Technology) as described in Knecht et al. (2017). Briefly, the LPR assay comprises 3-min alternating light and dark periods, for a total of four light:dark transitions. The first transition is counted as an acclimation period. The level was set at 525 LUX for the light period and integration time set to 6 s. Larval movement was recorded over each light:dark cycle, and average of total distance traveled per integration time point and area under the curve calculated.
Statistical analysis. For each concentration, n = 32 across two 96-well plates (16 animals per plate). The overall area under the curve for the final three light:dark transitions was computed for animals exposed to 0, 0.0625, or 0:125 ng=mL, and was compared with the 0 ng=mL TCDD concentration using a Kolmogorov-Smirnov test (p ≤ 0:01).
Statistical analysis. Each biological sample consisted of RNA from 20 pooled 48-hpf zebrafish with 4 biological replicates per condition (n = 4). Differential expression analysis followed the maintained Bioconductor workflow developed by , http://bioconductor.org/packages/release/workflows/vignettes/ RnaSeqGeneEdgeRQL/inst/doc/edgeRQL.html, and was conducted using R (version 3.5.1; R Development Core Team), RStudio (version 1.1.453; RStudio Team), and a custom R script provided as a supplementary file (Differential_expression_custom_script.R; Supplemental Material, S2. Differential expression custom R script; . The Bioconductor package edgeR (version 3.22.3) was used to normalize counts and identify differentially expressed genes McCarthy et al. 2012;Smyth 2007, 2008;Robinson and Oshlack 2010). Briefly, genes were filtered to exclude those with low counts across libraries, only keeping genes expressed in a minimum of four samples with average counts per million reads per sample above 0.82, which corresponds to a minimum read count of 10-20 (Y Lun et al. 2016). Filtered genes were then normalized across samples using the trimmed mean of M values (TMM) method to minimize composition bias between libraries (Robinson and Oshlack 2010). Differential expression of control vs. TCDD-exposed control or slincR morphants was determined with functions from edgeR, which uses the negative binominal generalized linear model extended by quasi-likelihood methods to fit the count data, the Cox-Reid profile-adjusted likelihood method to calculate dispersions, and empirical Bayes quasi-likelihood F-tests to calculate differential expression (Chen et al. 2014;Lun et al. 2016). The "robust = TRUE" option was used to protect the empirical Bayes estimates against the possibility of outlier genes with wideranging individual dispersions. Genes with a Benjamini-Hochberg (BH) adjusted p ≤ 0:05 were considered significantly differentially expressed. The biomaRt package (version 2.36.1) was used to connect Ensembl gene ID information to Ensembl BioMart annotation information (e.g., gene symbols, biotypes, human orthologs, etc.), and heatmaps were created using the R package ComplexHeatmap (version 1.18.1; Durinck et al. 2005;Durinck et al. 2009;Gu et al. 2016). Heatmap clustering was derived from TMM-normalized, regular log-transformed gene values scaled by z-score (Love et al. 2016). To understand the functional consequences of TCDD exposure and slincR knockdown, we performed biological process network enrichment analysis and GO term enrichment on the significant differentially expressed human orthologs using GeneGo MetaCore (version 6.31 build 68930) from Clarivate Analytics as described in Haggard et al. (2016). Only enriched biological process networks and GO terms with a false discovery rate (FDR) adjusted p ≤ 0:05 were considered significant. Sequencing data and processing details have been deposited in the NCBI Gene Expression Omnibus (GSE106131).
Cartilage staining and lower jaw cartilage morphometrics. Zebrafish embryos were microinjected with control and slincR MOs and exposed to 0.1% DMSO or 1 ng=mL TCDD, as described above. At 72 hpf, larvae were euthanized by overdose of buffered tricaine methanesulfonate (MS-222, 200-300 mg=L) by prolonged immersion and monitored under a dissecting microscope until the hearts stopped beating (approximately 15 min). The animals were fixed in 4% paraformaldehyde overnight at 4°C. Pigmentation was removed by incubating fixed larval samples for 1 h in a mixture of 3% H 2 O 2 =1% KOH. The cartilage was stained with 0.4% Alcian Blue 8GX (Sigma-Aldrich) in 70% ethanol and 80 mM MgCl 2 as described in Walker and Kimmel (2007). Briefly, fixed animals were washed with PBS and dehydrated with an ethanol gradient. After removal of ethanol, the cartilage stain was added to larvae and rocked overnight at room temperature. The samples were washed with water and cleared using glycerol and KOH solutions. Larvae were imaged in 0.8% low-melt agarose at room temperature using a Keyence BZ-X700 at 10X with 0.45 aperture and processed with the BZ-X Analyzer software. Morphometric analysis (n = 9-10 embryos) of ventral larval pharyngeal cartilages was performed using ImageJ v1.51j8 and two customized macros (Customized_Set_Origin_Tool.ijm; Customized_Click_Coordinates_Tool.ijm; Supplemental Material, S3. Customized Set Origin Tool and S4. Customized Click Coordinates Tool), as described in Xiong et al. (2008).
Hemorrhage screen. Zebrafish embryos were microinjected with control and slincR MOs and exposed to 0.1% DMSO or 1 ng=mL TCDD, as described above. At 48 hpf, the embryos were evaluated under a dissecting microscope for the presence or absence of hemorrhaging.
Statistical analysis. The DMSO and TCDD samples consisted of 3 or 5 biological replicates, respectively. Each biological replicate contained 10-12 individual 48-hpf zebrafish. Results were statistically analyzed with R (version 3.4.1; R Development Core Team) in the RStudio IDE (version 1.0.143; RStudio Team). The data were initially assessed for equal variance and normality using Levene's test and the Shapiro-Wilk test, respectively. Statistical significance (n = 3 or 5, p < 0:001) was determined using the Type III one-way ANOVA for unbalanced experimental designs from the car package (version 3.0-0) with a Tukey post hoc test from the multcomp package (version 1.4-8), and the data were graphed using ggplot2 (version 3.0.0; Fox and Weisberg 2011;Hothorn et al. 2008;Wickham 2016). The experiments were independently repeated a minimum of two times.
PART III: Investigation of PAHs That Increase slincR Expression 16 PAH RNA-seq sample preparation and analysis. The zebrafish embryos were exposed as described in Geier et al. (2018). Briefly, the chorions were enzymatically removed from 4-hpf zebrafish using a custom automated dechorionator. At 6 hpf, the embryos were placed into a 96-well round-bottom plate using automated embryo placement robots (Mandrell et al. 2012). Each well contained a single embryo and 100 lL embryo media. The 16 PAHs (Table 1) were analytical grade standards obtained from AccuStandard, Chiron Chemicals, and Santa Cruz Biotechnology. To make stock solutions, the 16 PAHs were dissolved in 100% DMSO. See Geier et al. (2018) for detailed stock information. A Hewlett Packard D300e chemical dispenser was used to expose the embryos in 96-well plates. To achieve optimal solution uniformity, the plates were sealed with Parafilm ® , covered in foil, and shaken overnight at 235 rpm in an orbital shaker at 28°C (Truong et al. 2016). The exposure plates were then transferred to an incubator at 28°C for the duration of the exposure. The EC 80 's of the 16 PAHs were experimentally determined by the individual compound response. The 16 PAHs were first assessed at a wide range of concentrations (50, 35.6, 11.2, 5, and 1 lM, n = 32). Of the 16 PAHs, 9 did not result in morphological malformations and were subsequently tested at only 50 lM. The remaining 7 were tested using a definitive concentration range (11 concentrations, n = 24) determined between the highest concentration to not elicit any morphological effect and the lowest concentration that resulted in near 100% morbidity or mortality. The final EC 80 concentrations tested are listed behind each chemical name in Table 1. Sigmoidal curves are often fit to concentration-response data. Here, we used a Hill Model (specifically a four-parameter loglogistic function) that was fit to the mean percentage of affected individuals for any morphological end point measured. All curves were fit with the drm() function from the drc package in R. This function uses least squares estimation to fit the curves. The Hill model was applied to estimate a concentration that caused 80% effects (EC 80 ). The computed EC 80 was confirmed prior to conducting the RNA-seq. The data for the morphological effects of the 16 PAHs are reported in Geier et al. (2018). Total RNA was isolated and quantified from pooled groups of eight 48-hpf zebrafish as described in the "RNA Extraction" section of the "Methods" section. Each treatment group contained 4 biological replicates. Total RNA samples were sent to Oregon State University Center for Genome Research and Biocomputing Core facilities for library preparation and sequencing. The 100-bp single-end read libraries were prepared using the Wafergen Robotic PolyA Enrichment Library Prep and Wafergen Robotic Stranded RNA Library Prep Kits. The samples were randomized across 6 lanes and sequenced using the Illumina HiSeq ® 3000.
Statistical analysis. The reads were processed and analyzed as described in the "TCDD RNA-seq Sample Preparation and Analysis section" of the "Methods" section. Each biological sample consisted of RNA from 8 pooled 48-hpf zebrafish with 4 biological replicates per condition (n = 4, BH-adjusted p ≤ 0:05).

PART IV: Identification of Potential slincR Mouse and Human Orthologs
Mouse TCDD-exposed dataset. All procedures were approved by the University of Wisconsin Animal Care and Use Committee and conducted in accordance with the National Institutes of Health (NIH) Guide for Care and Use of Laboratory Animals. Pregnant C57BL/6J mice were treated with 5 lg=kg TCDD on E13.5 following previously established protocols (Lin et al. 2003). At E16.75, dams were euthanized via CO 2 asphyxiation, and the urogenital sinuses (UGS) were dissected from fetuses as previously described (Branam et al. 2013). Each UGS was immediately placed into a 1:5-mL microfuge tube containing 300 lL of 1% trypsin (Difco, 215240) in PBS and incubated on ice for 30 min. Collagenase (Sigma C9891) was added to a final concentration of 1 mg=mL, followed by an additional 30-45 min incubation on ice. A dissecting microscope was used to mechanically separate UGS from fetal urogenital epitheliums (UGEs) after which the bladder and urethral epithelium were removed, leaving only UGE. Each treatment consisted of 3-4 biological replicates consisting of either 5 or 6 individual UGEs that were pooled for RNA isolation. Total RNA was purified from each UGE using the RNeasy ® system (Qiagen) and analyzed using the Agilent Bioanalyzer 2100 and an Agilent RNA 6000 Pico Kit (Agilent Technologies). The samples were sequenced at the University of Wisconsin Madison Biotechnology Center, using an Illumina HiSeq ® 2500. Processing of NGS was performed as described above in the "RNA sequencing analysis"' section with the following modifications: The genome index files were generated with the GRCm38 (Ensembl Release 87) genome file, ftp://ftp.ensembl.org/ pub/release-87/fasta/mus_musculus/dna/Mus_musculus.GRCm38. dna.toplevel.fa.gz; reads were trimmed using the FASTX-Toolkit (version 0.0.13; http://hannonlab.cshl.edu/fastx_toolkit/index.html) with options: fastx_trimmer -t 50 -Q 33; TopHat2 options were -library-type fr-unstranded -no-mixednum-threads 10; and gene counts were determined using the htseq-count script from HTSeq with the Ensembl Release 87 GTF annotation, ftp://ftp.ensembl.org/pub/release-87/gtf/mus_ musculus/Mus_musculus.GRCm38.87.gtf.gz and options: -f bam -i gene_id -m intersection-nonempty -s no.
Statistical analysis. Each biological sample consisted of RNA from 5-6 pooled UGE tissue with 3-4 biological replicates per condition. The differential expression analysis was performed as described above in the RNA sequencing analysis section (n = 3-4, BH-adjusted p ≤ 0:05).
Identification of potential mammalian slincR orthologs. The slncky Evolution Browser contains alignments and evolutionary metrics of lncRNAs conserved in the mouse (mm10, GRCm38) and human (hg38, GRCh38) genomes and was used to identify the potential human ortholog of slincR . We performed pairwise sequence alignments using the online global alignment tool Needle (EMBOSS) using the default settings (https://www.ebi.ac.uk/Tools/psa/). To compare tissue-specific expression, we downloaded RNA sequencing data from the National Center for Biotechnology Information BioProject database for the mouse, 2610035D17Rik (BioProject: PRJNA66167; Gene ID: 72386), and human, LINC00673 (BioProject PRJEB4337; Gene ID: 100499467), orthologs of slincR, as well as for sox9

PART I: Transcriptional Regulation by slincR
Use of qRT-PCR CHART to evaluate slincR enrichment at the sox9b locus. The slincR probe was shown to be capable of enriching the slincR transcript ( Figure 1A) and was shown to be specific as demonstrated by the absence of the slincR transcript band in the sense probe samples of a representative 1.2% agarose gel ( Figure 1B). We tested multiple regions of the sox9b promoter, and only the 5 0 UTR of the promoter was consistently enriched in all three biological replicates in both the DMSO-and TCDD-treated samples ( Figure 1C). The 5' UTR is also present in the sox9b transcript, so we tested the same region in CHARTisolated RNA; however, the RNA samples did not result in detectable amplification (no measurable enrichment).
To determine if sox9b expression levels have an effect on the expression of slincR, we microinjected single-cell zebrafish embryos with a splice-blocking morpholino targeting sox9b (sox9b-MO) or with a control morpholino (Con-MO) and collected RNA at 48 hpf. We used qRT-PCR to measure the relative expression levels of sox9b and slincR. The sox9b-MO reduced the levels of sox9b transcripts by 71% and resulted in a significant increase in slincR expression when compared with control morphants ( Figure 1D).
To determine if a concentration-response relationship exists between the gene expression of sox9b and slincR in response to TCDD exposure, we developmentally exposed embryonic zebrafish to 0, 0.0625, 0.125, 0.25, 0.5, or 1:0 ng=mL of TCDD and measured the relative expression levels of sox9b and slincR in 48-hpf wild-type zebrafish. We also measured cyp1a expression as a marker of AHR activation. The expression levels of slincR and cyp1a were significantly elevated beginning at the 0:0625 ng=mL TCDD concentration and increased with increasing concentrations of TCDD, in comparison with 0 ng=mL TCDD (i.e., 0.1% DMSO; Figure 2A). The expression levels of sox9b were significantly decreased in the 0.5 and 1:0 ng=mL TCDD concentrations, in comparison with 0 ng=mL TCDD (Figure 2A).
In addition to gene expression changes, we also examined the embryos for developmental toxicity when exposed to the six concentrations of TCDD (0-1 ng=mL) using a 17 end point morphology screen at 120 hpf, as described in Truong et al. (2014). Developmental exposures of 0.25, 0.5, or 1:0 ng=mL TCDD induced significant malformations across multiple end points, including edemas (of the yolk and pericardium) and jaw, snout, eye, and pectoral fin malformations ( Figure 2B, C). Developmental exposures of 0.0625 or 0:125 ng=mL did not induce a significant level of physical malformations ( Figure 2B, C); however, both concentrations resulted in a significant hyperactive phenotype in a larval photomotor response (LPR) assay in comparison with 0 ng=mL TCDD ( Figure  2D). Larval zebrafish exhibiting any visible malformation were necessarily excluded, which prevented the analysis of the 0.25, 0.5, or 1:0 ng=mL TCDD concentrations in the LPR assay.
PART II: slincR Contributions to TCDD Toxicity Unbiased transcriptome profiling in TCDD-exposed zebrafish morphants. To gain insight into the functional role of slincR in TCDD-induced toxicity, we performed whole embryo transcriptome profiling on 48-hpf slincR and control morphants exposed to 0.1% DMSO or 1 ng=mL TCDD. In control morphants exposed to TCDD, we identified 132 genes that were differentially expressed when compared with vehicle control (n = 4, BH-adjusted p ≤ 0:05, Figure 3A and Excel Table S1). In slincR morphants exposed to TCDD, we identified 80 differentially expressed genes when compared with vehicle control (n = 4, BH-adjusted p ≤ 0:05, Figure 3A and Excel Table S2). Over 25% (22) of these genes were unique to slincR morphants, whereas 72.5% (58) overlapped with the control morphant differentially expressed gene list ( Figure 3A). In the slincR and control morphant differentially expressed gene lists, 90% (72) and 81% (107) of the genes had an increase in expression in response to TCDD exposure, respectively. As expected, the sox9b transcript was only significantly decreased in the control morphant (TCDD-DMSO) dataset (Excel Table S1 and S2).
To visualize the expression profile, we performed bidirectional hierarchical clustering of the differentially expressed genes . Each condition had 3 biological replicates, where 1 replicate consisted of approximately 500 48-hpf zebrafish. The qRT-PCR data were first normalized to 1% input control, such that 6.644 cycles (i.e., dilution factor log 2 ð100Þ) was subtracted from the cycle threshold value (CT; i.e., number of PCR replication cycles required for the sample signal to exceed background levels) of the diluted input and used to calculate the DCT for the two probe sets (DCT = CT½probe − CT½1% input-6:644). The RNA yield was calculated using the following equation (2 −DDCT × 100%). We assigned samples that did not amplify (no enrichment) a CT value of 40. (B) Representative slincR qRT-PCR CHART products from panel (A) run on a 1.2% agarose gel. (C) qRT-PCR CHART enrichment of slincR RNA at multiple positions (−2042 bp, −963 bp, and −502 bp) downstream of the sox9b transcription start site and 5' untranslated region. Each condition had 3 biological replicates (n = 3), where 1 replicate consisted of approximately 500 48-hpf zebrafish. Expression values were normalized to 1% input control as described for panel (A), except for DNA fold enrichment. We next adjusted relative to the sense probe (DDCT = DCT½slincR-probe − DCT½sense-probe), and then fold enrichment was calculated (2 −DDCT ). We assigned samples that did not amplify (no enrichment) a CT value of 40. (D) qRT-PCR relative expression of slincR and sox9b mRNA in 48-hpf whole embryo control and sox9b morphants. Expression values were analyzed with the 2 −DDCT method and normalized to b-actin, whereas the control morphants served as the calibrator. Each sample represents a pool of 20 embryos, each condition had a minimum of four biological replicates (n = 4), and the data were analyzed using a one-way ANOVA with a Tukey post hoc test (p < 0:01 in comparison with control morphant = Ã Ã). Error bars (A, C, and D) indicate standard error of the mean. Note: CHART, capture hybridization analysis of RNA targets; DMSO, dimethyl sulfoxide; hpf, hours post fertilization; TCDD, 2,3,7,8-tetrachlorodibenzo-p-dioxin; UTR, untranslated region.   Relative Expression Figure 2. Analysis of the concentration-response effects of developmental exposure (1-h exposure at 6 hpf) to TCDD on gene expression at 48 hpf and morphological malformations at 120 hpf. (A) qRT-PCR relative expression of cyp1a, slincR, and sox9b transcripts in 48-hpf wild-type embryos exposed to 0, 0.625, 0.125, 0.25, 0.5, or 1:0 ng=mL TCDD. For all assays, 0.1% DMSO served as the vehicle control and is listed as 0 ng=mL TCDD. Expression values were analyzed with the 2 −DDCT method and normalized to b-actin using the 0 ng=mL TCDD concentration as the calibrator. Each experimental unit represents a pool of 20 embryos, and each treatment group included four biological replicates (n = 4). The data were log 2 -transformed and analyzed using a one-way ANOVA with a Dunnett post hoc test (p < 0:001 in comparison with 0 ng=mL TCDD = Ã ÃÃ). All error bars indicate standard error of the mean. (B) Evaluation of 17 physical malformations at 120-hpf on wild-type zebrafish exposed to six concentrations of TCDD (0-1 ng=mL) across two 96-well plates. Nonsignificant malformations (otic vesicle, somite, circulation, pigmentation, swim bladder, notochord distortion, and alterations in touch response) are excluded. The horizontal axis displays the TCDD concentrations tested, and the malformation examined is listed above each box. The incidence across all replicates is plotted as stacked points. For each malformation, the stacked points exceeding the binomial significance threshold are represented in light gray (top stack). The data were analyzed using a Fisher's exact test with a Bonferroni correction for multiple comparisons (n = 32, p < 0:01). (C) Representative lateral images of 120-hpf zebrafish for each concentration of TCDD tested. The bar in top left corner indicates 100 lM. (D) Larval photomotor response (LPR) in 120-hpf wild-type embryos developmentally exposed to 0, 0.0625, or 0:125 ng=mL TCDD equally divided across two 96-well plates using the ViewPoint ZebraBox larvae screening system. For each concentration of TCDD, the overall area under the curve was analyzed for the last 3 light:dark cycles in comparison with 0 ng=mL TCDD using a Kolmogorov-Smirnov test (n = 32, p ≤ 0:01 with the top 30 largest log 2 ðfold changesÞ in control ( Figure 3B) and slincR ( Figure 3C) morphants. TCDD exposure produced a robust transcriptional profile as indicated by the two primary transcript clusters that separate based on treatment status (blue and red bars at the top of the heatmaps).
To relate the observed transcriptional changes to human health and understand the functional consequences of TCDD exposure and slincR knockdown, we identified orthologous human Ensembl genes and performed biological process network enrichment and GO term enrichment analyses on the converted gene lists using MetaCore GeneGo software (FDR-adjusted p ≤ 0:05) as described in Haggard et al. (2016). The MetaCore process network enrichment for the slincR morphant dataset was limited because of its small size but included four networks involved in immune response, regulation of epithelial-to-mesenchymal transition, ESR1-nuclear pathway, and reproduction, that were also  significantly enriched in the control morphant data (Table 2 and 3). We observed multiple enriched process networks involved in cartilage, bone, and blood vessel development in the control morphant dataset. The GO term enrichment analyses resulted in about half of GO-enriched terms overlapping between Con-MO and slincR-MO datasets (Excel Table S3 and S4). We noticed that several cartilage and skeletal developmental processes were uniquely enriched in the control morphant data (Table 4), whereas angiogenesis and vasculature developmental processes were uniquely enriched in the slincR dataset (Table 5).
Effect of slincR knockdown on jaw development after treatment with TCDD. To investigate the phenotypic impact of slincR expression on TCDD-induced jaw malformations, we developmentally exposed slincR and control morphants to 0.1% DMSO or 1 ng=mL TCDD for 1 h at the shield stage ( ∼ 6 hpf) and evaluated the cartilage of 72-hpf zebrafish larvae. Both the control and slincR morphants exposed to TCDD had a 100% incidence of cartilage malformations at 72 hpf. The cartilage of DMSOexposed slincR and control morphants was indistinguishable; however, slincR morphants exposed to TCDD had a significant difference in the relative position of landmark structures in comparison with TCDD-exposed control morphants (n = 9 or 10, p < 0:05; Figure 4A, B). Unexpectedly, knocking down slincR expression in TCDD-exposed (Ahr2-activated) zebrafish embryos resulted in an abnormal junction between hyosymplectic and ceratohyal cartilages in comparison with control morphants exposed to TCDD ( Figure 4B, region of interest indicated by arrows in trace column).
Effect of slincR knockdown on TCDD-induced hemorrhaging. To elucidate the phenotypic impact of slincR expression on TCDD-induced hemorrhaging, we developmentally exposed slincR and control morphants to 0.1% DMSO (vehicle control) or 1 ng=mL TCDD for 1 h at the shield stage ( ∼ 6 hpf). At 48 hpf, the embryos were evaluated for the presence or absence of hemorrhaging as shown in (Figure 5A, B; blood pooling indicated by arrows). In control and slincR morphants, exposure to DMSO did not result in hemorrhaging ( Figure 5A). In both the control and slincR morphants, exposure to TCDD resulted in a significant increase in the percent incidence of zebrafish displaying a hemorrhaging phenotype in comparison with their respective vehicle controls (n = 3 or 5, p < 0:001; Figure 5A); however, the slincR morphants exposed to TCDD had a significant decrease in the percentage of animals with the hemorrhaging phenotype in comparison with control morphants exposed to TCDD ( Figure 5A).

PART III: Investigation of PAHs That Increase slincR Expression
slincR is up-regulated by multiple environmentally relevant PAHs. To screen for environmentally relevant PAHs that activate the AHR signaling pathway and up-regulate slincR expression, we mined an unpublished 16 PAH RNA-seq dataset from 48-hpf whole zebrafish embryos generated by Dr. M. Geier from the Tanguay lab. Of the 16 PAHs screened, six were associated with a significant increase in cyp1a and slincR expression, of which three were from the U.S. Environmental Protection Agency's priority PAH list (Table 1). Three additional PAHs were associated with a significant increase in cyp1a expression alone, of which two were from the EPA's priority PAH list. None of the 16 PAHs were associated with a significant decrease in sox9b expression.

PART IV: Identification of Potential slincR Mouse and Human Orthologs
To identify the potential mouse ortholog of slincR, we mined an unpublished RNA-seq dataset from male and female mouse E16.75 urogenital epithelial tissue exposed to 5 lg=kg of TCDD, Table 3. Significantly enriched MetaCore process networks from 48-hpf slincR-MO zebrafish exposed to 1 ng=mL TCDD for 1 hour at the shield stage, in comparison with the vehicle control (0.1% DMSO).  Note: The Bioconductor package edgeR (statistical methodology based on the negative binomial distribution) was used to calculate the significant differentially expressed genes (n = 4, BH-adjusted p ≤ 0:05), in comparison with the vehicle control (0.1% DMSO). The human orthologs of the gene list were then submitted to MetaCore to identify significantly enriched process networks using a hypergeometric distribution, where the p-value is the probability that a gene set maps to a manually curated GeneGo Process Network or is overrepresented in comparison with the background gene list. Enriched process networks were considered significant with an FDR ≤ 0:05. BH, Benjamini-Hochberg; BMP, bone morphogenic protein; Con-MO, control morphant; DMSO, dimethyl sulfoxide; ECM, extracellular matrix; EMT, epithelial to mesenchyme transition; FDR, false discovery rate adjusted p-value; FSH-beta, follicle-stimulating hormone beta subunit; GDF, growth differentiation factor; hpf, hours post fertilization; slincR-MO, slincR morphant; TCDD, 2,3,7,8-tetrachlorodibenzo-p-dioxin.
which was generated by the Richard Peterson lab from the University of Wisconsin-Madison. In the zebrafish genome, slincR is located adjacent and antisense to sox9b; therefore, we searched for a lncRNA that was significantly increased in response to TCDD and had a conserved genomic position and orientation relative to Sox9 in the mouse genome. We identified a single lncRNA (2610035D17Rik) that matched these criteria ( Figure 6A, B). We used the slncky Evolution Browser to search for the human ortholog of mouse 2610035D17Rik and identified LINC00673 ( Figure  6C; . A sequence comparison of the zebrafish slincR sequence (466 bp) with the mouse (1802 bp) and human (2275 bp) lncRNAs resulted in a 17.6% and 13.4% sequence identity, respectively. The mouse and human lncRNAs shared 43.1% sequence identity. To determine if the mouse and human lncRNAs have similar tissue-specific expression patterns in comparison with slincR in zebrafish, we downloaded mouse and human RNA-seq expression data from NCBI and graphed the expression of the potential slincR orthologs and Sox9/SOX9 across multiple tissues ( Figure 6D, E). Concordant with slincR, expression of both 2610035D17Rik and LINC00673 was found in the brain and central nervous system, and 2610035D17Rik expression was also identified in the developing limb bud.

Discussion
The AHR is required for proper vertebrate development and homeostasis; however, activation of the receptor by ubiquitous environmental pollutants, such as PAHs, can lead to adverse developmental and cognitive effects in humans and wildlife (Carney et al. 2006;Schneider et al. 2014). Dysregulation of the AHR signal transduction pathway is associated with multiple diseases, including prostate and coronary artery disease (Huang et al. 2015;Schneider et al. 2014;Vezina et al. 2009). AHR dysregulation has also been implicated in many cancer types and mediates a number of steps in tumor progression (Murray et al. 2014;Opitz et al. 2011). Recent publications suggest the AHR may play a role in maintaining the cancer stemlike phenotype and may exert transcription-independent functions to mediate resistance to treatment in adenocarcinoma and nonsmall cell lung cancers, respectively (Yan et al. 2018;Ye et al. 2018). A better understanding of the molecular mechanisms that lead to toxic outcomes will facilitate the necessary shift from observations of apical end points (e.g., malformations, mortality) to predicting the potential of a chemical to alter or interfere with biologically conserved pathways . We Table 4. Significantly enriched unique MetaCore GO processes related to skeletal and cartilage development from 48-hpf Con-MO zebrafish exposed to 1 ng=mL TCDD for 1 h at the shield stage, in comparison with the vehicle control (0.1% DMSO).
Biological process FDR Associated human gene Note: The Bioconductor package edgeR (statistical methodology based on the negative binomial distribution) was used to calculate the significant differentially expressed genes (n = 4, BH-adjusted p ≤ 0:05), in comparison with the vehicle control (0.1% DMSO). The human orthologs of the gene list were then submitted to MetaCore to identify significant biological processes using a hypergeometric distribution, where the p-value is the probability that a gene set maps to a manually curated GO process or is overrepresented in comparison with the background gene list. Enriched GO processes were considered significant with an FDR ≤ 0:05. BH, Benjamini-Hochberg; Con-MO, control morphant; DMSO, dimethyl sulfoxide; FDR, false discovery rate adjusted p-value; GO, gene ontology; hpf, hours post fertilization; TCDD, 2,3,7,8-tetrachlorodibenzo-p-dioxin. Table 5. Significantly-enriched unique MetaCore GO processes related to angiogenesis and vasculature development from 48-hpf slincR-MO zebrafish exposed to 1 ng=mL TCDD for 1 h at the shield stage, in comparison with the vehicle control (0.1% DMSO).
previously used the zebrafish model to identify an ahr2-dependent long noncoding RNA, slincR, whose expression increased upon exposure to multiple Ahr2 ligands . We demonstrated that slincR was required for the proper expression of sox9b during development and sox9b repression in response to TCDD, a strong AHR ligand ).

PART I: Transcriptional Regulation by slincR
Our data suggested that slincR acted in cis to repress sox9b expression. Biochemical evidence has demonstrated that many lncRNAs act as guides to chromatin-modifying enzymes and/or platforms for protein complexes to regulate the three-dimensional structure of the genome (Quinn and Chang 2016). In the present study, our data suggest that slincR regulates the expression of sox9b via enrichment at the 5' UTR of the sox9b locus. Further experiments are needed to distinguish whether slincR binds directly to the DNA or indirectly via a protein/RNA intermediate. Surprisingly, the DMSO-and TCDD-exposed samples produced similar yields of slincR RNA and sox9b 5' UTR DNA enrichment, despite the significant increase in slincR expression induced upon exposure to TCDD. One possible explanation is that the probe accessibility of the slincR transcript is altered due to changes in the interacting molecular partners, which may be condition specific. Additionally, knocking down sox9b expression resulted in a significant increase in slincR expression. Our previous publication demonstrated that slincR was required for the TCDD-induced repression of sox9b, had a significant effect on the expression of known sox9b target genes, and is expressed in tissues with sox9b essential functions . These data imply that sox9b and slincR may share overlapping regulatory networks. We detected a clear concentration-response relationship between TCDD exposure and the gene expression of slincR and sox9b. We also evaluated the morphology and behavior of exposed zebrafish in order to correlate the gene expression with the observed toxicity. Although the 0.25, 0.5, and 1:0 ng=mL concentrations of TCDD produced visible malformations in the eye, brain, jaw, snout, and other tissues, only the two highest concentrations of TCDD tested (0.5 and 1:0 ng=mL) were able to both significantly increase slincR and significantly decrease sox9b expression levels. In response to developmental exposure to 1 ng=mL TCDD, sox9b expression is reduced in the eye, brain, otic vesicle, jaw region, and snout region . For the 0:25 ng=mL TCDD concentration, the discrepancy . TCDD-induced hemorrhaging in control and slincR morphants. Control (Con-MO) and slincR morphants (slincR-MO) were developmentally exposed to 0.1% DMSO or 1 ng=mL TCDD, and 48-hpf zebrafish embryos were evaluated for the presence or absence of hemorrhaging. (A) The percent incidence of the hemorrhaging phenotype in control and slincR morphants exposed to DMSO or TCDD. The DMSO and TCDD samples consisted of 3 and 5 biological replicates, respectively. Each biological replicate contained 10-12 individual 48-hpf zebrafish. All error bars indicate standard error of the mean. Statistical significance was determined using the Type III one-way ANOVA for unbalanced experimental designs from the car package with a Tukey post hoc test (n = 3 or 5, p < 0:001 = Ã ÃÃ). (B) Representative images of 48-hpf wild-type zebrafish exposed to DMSO (normal phenotype) or TCDD (hemorrhaging phenotype). The bar in the top right corner indicates 100 lM and arrows point to sites of hemorrhaging. Note: DMSO, dimethyl sulfoxide; hpf, hours post fertilization; TCDD,2,3,7, between the nonsignificant sox9b gene expression and significant malformations in the cartilage area may be because the RNA was isolated from a pool of whole embryos, which does not allow for tissue-specific resolution.
PART II: slincR Contributions to TCDD Toxicity SOX9 is required for proper vertebrate development and regulates cell maintenance and specification during adulthood (Barrionuevo et al. 2016;Furuyama et al. 2011;Jo et al. 2014). Overexpression of SOX9 is associated with liver fibrosis and multiple cancer types (Pritchett et al. 2011). In zebrafish, sox9b has a causal role in the Ahr2 toxicity pathway. Antisense knockdown of sox9b is sufficient to recapitulate the craniofacial cartilage malformations induced by exposure to TCDD (Xiong et al. 2008). Developmental TCDD-induced repression of sox9b decreases both the size and number of chondrocytes to produce malformed craniofacial cartilages (Burns et al. 2015). In addition, ahr2-mutant zebrafish also display craniofacial skeletal abnormalities, further supporting the intersection of the AHR and SOX9 regulatory networks (Garcia et al. 2018;Goodale et al. 2012). To understand the functional role of slincR in TCDD-induced developmental toxicity, we investigated global transcriptional responses and performed functional enrichment analyses on control and slincR morphants   Figure 6. Identification of the potential slincR mouse and human orthologs. We identified the potential mouse ortholog of slincR (2610035D17Rik) based on (A) the conserved genomic location and orientation relative to Sox9 and (B) the significant increase in expression in TCDD-exposed mouse urogenital tissue samples from embryonic day 16.5 (F, female and M, male). Each biological sample consisted of RNA from 5-6 pooled mouse urogenital tissue with 3-4 biological replicates per condition. Error bars indicate standard error of the mean. The image from (A) was downloaded from the Ensembl mouse genome (GRCm38), and Sox9 and 2610035D17Rik are highlighted. (C) Downloaded and formatted image from the slncky Evolution Browser of the mouse (2610035D17Rik) and human (LINC00673) conserved lncRNA orthologs. To determine the tissue-specific expression of (D) mouse-(2610035D17Rik) and (E) human-conserved (LINC00673) lncRNAs relative to Sox9/SOX9, we downloaded RNA-seq expression data from NCBI BioProjects PRJNA66167 and PRJEB4337, respectively. Error bars indicate standard error of the mean. Note: CPM, counts per million; RPKM, reads per kilobase million. exposed to 1 ng=mL TCDD. Upon exposure to TCDD, an increase in slincR expression is required for the functional enrichment of many sox9b-regulated processes, including cartilage development, ossification, and bone remodeling, as well as cell adhesion and cell matrix interactions (SOX9 functions reviewed in Jo et al. 2014). We saw evidence that slincR knockdown in the presence of TCDD can also alter the structure of craniofacial cartilage. In our morphometric analysis of the cartilage structure, both the control and slincR morphants exposed to TCDD exhibited cartilage malformations; however, the slincR morphants showed an abnormal junction between hyosymplectic and ceratohyal cartilages in comparison with control morphants. The functional enrichment and morphometric analyses support the hypothesis that, upon TCDD exposure, Ahr2 increases the expression of slincR, which in turn represses sox9b via enrichment at the 5' UTR, to produce craniofacial cartilage malformations. It is possible that the abnormal cartilage morphology in the slincR morphants is due to incomplete slincR knockdown and thus incomplete sox9b repression (i.e., higher expression of sox9b) relative to the levels in the control morphant samples exposed to TCDD. In other words, differing expression levels of sox9b may be responsible for the altered cartilage structure, because sox9b is a master regulator of cartilage development.
In addition to the AHR-dependent induction of slincR, other factors may be acting to repress sox9b expression. For example, crosstalk between AHR and Wnt=b-catenin signaling pathways is well established (Schneider et al. 2014;Wincent et al. 2015). In mammalian epithelial stem cell in vitro experiments, SOX9 and WNT signaling cooperate in a mutually repressive manner to regulate proliferation, differentiation, and quiescence (Menzel-Severing et al. 2018). In mice, deletion of b-catenin in head mesenchyme provided in vivo evidence that suggests b-catenin is essential for proper skeletal linage differentiation by inhibiting mesenchymal osteoblastic cells from entering the chondrocyte linage (Hill et al. 2005). Gene expression measurements in the b-catenin mutant displayed a significant increase in expression of Sox9 and its upstream regulator Runx2, which argues that b-catenin is not directly regulating Sox9. In zebrafish, several publications have supported the hypothesis that xenobiotic activation of the AHR leads to a disruption in Wnt=b-catenin signaling, and that this disruption in Wnt signaling has a causal role in AHR-mediated toxicities (Mathew et al. 2008;Wincent et al. 2015). Using the zebrafish caudal fin regeneration model, we have previously shown that TCDD-induced AHR activation blocks the tissue regeneration process via activation of the Wnt=b-catenin signaling pathway, leading to improper formation of the wound epithelium and blastema (Mathew et al. 2008). The TCDD-induced block in caudal fin regeneration is prevented by antisense knockdown of a Wnt coreceptor (lrp6) that induces b-catenin signaling. In addition, exposure to a pharmacological agent (6-bromoindirubin-3 0 -oxime) that overactivates the Wnt signaling pathway phenocopied the TCDD-induced block in fin regeneration. In response to AHR activation, it is likely that genes involved in multiple signaling pathways are disrupted and acting in concert to repress expression of sox9b and/or inhibit normal cartilage development.
The AHR regulates vasculature remodeling of the developing embryo, and dysregulation can lead to abnormal development of vasculature structures (Lahvis et al. 2000;Walisser et al. 2004). Developmental exposure to TCDD or PAHs can produce blood circulation defects and cyp1a1 expression in the vasculature (Andreasen et al. 2002;Chlebowski et al. 2017;Henry et al. 1997). Our transcriptome profiling and functional enrichment analyses suggest slincR expression regulates processes involved in angiogenesis, vascular smooth-muscle cell proliferation, blood vessel endothelial cell migration, blood coagulation, and blood vessel morphogenesis.
Upon exposure to a strong Ahr2 ligand, a toxicity phenotype that occurs early in zebrafish development is disruption of blood cell development and hemorrhaging (Belair et al. 2001;Henry et al. 1997). To elucidate the phenotypic impact of slincR expression on TCDD-induced hemorrhaging, we developmentally exposed control and slincR morphants to TCDD and recorded the presence or absence of hemorrhaging. In concordance with the RNA-seq data, we demonstrated that upon TCDD exposure, preventing the induction of slincR can prevent a significant percentage of zebrafish from displaying the hemorrhaging phenotype; however, knocking down slincR did not result in the complete absence of the hemorrhaging phenotype. Although this discrepancy could be due to the incomplete knockdown of slincR, it is also possible that the activation of the AHR results in the dysregulation of multiple regulatory networks that are required for normal angiogenesis and vascular development. Our data suggest ahr2dependent expression of slincR plays a causal role in TCDDinduced hemorrhaging. Additional experiments are required to understand how slincR expression regulates processes such as angiogenesis and vasculature development at the gene-network, cellular, and tissue levels.

PART III: Investigation of PAHs that Increase slincR Expression
The AHR is ligand-activated by numerous PAHs, which pose a major concern due to their potential toxic effects on ecosystems and human health (Chlebowski et al. 2017;Geier et al. 2018;Perera 1997). For example, exposure of pregnant women to high levels of air pollution, which is a complex heterogeneous mixture containing numerous PAHs, has been associated with an increase in the risk of impaired neural development and congenital heart defects in their children (Brook et al. 2004;Edwards et al. 2010). To screen for environmentally relevant PAHs that upregulate slincR expression and activate the AHR signaling pathway, we mined a 16 PAH RNA-seq dataset from 48-hpf whole zebrafish embryos generated by Dr. M. Geier from the Tanguay lab. We identified six PAHs that significantly up-regulate slincR expression, including benzo[k]fluoranthene, dibenzo[a,h]pyrene, and dibenzo[a,i]pyrene from the U.S. Environmental Protection Agency's priority PAH list. Our screen also identified three PAHs (retene, dibenzo[a,h]pyrene, and dibenzo[a,i]pyrene) that cause a significant increase in slincR expression and were recently shown to induce cyp1a1 vasculature expression (Geier et al. 2018). The Geier study also reported that fluoranthene and 9-methylanthracene induced cyp1a1 expression in the vasculature; however, our data did not show a significant increase in slincR expression for these two PAHs, emphasizing the ligand-specific transcriptional responses of AHR activation (Geier et al. 2018).
None of the PAHs screened resulted in a significant decrease in sox9b expression even though benzo[k]fluoranthene and benzo [j]fluoranthene were able to produce larger log 2 ðfold changeÞ increases (4.0 and 3.9) in slincR expression than did TCDD exposure (3.1). LncRNAs have a more restricted tissue-specific expression pattern in comparison with mRNAs (Fatica and Bozzoni 2014). Thus, it is possible that slincR induction was localized to tissue where sox9b is not highly expressed. We previously showed that exposure to TCDD increased slincR expression in the otic vesicle and jaw or snout regions . It is also possible that slincR regulates additional genes, beyond sox9b. Future experiments will identify the tissues in which slincR expression is increased upon exposure to select PAHs. Elucidating the transcriptional responses and target organs of ligand-specific AHR activation is a key process required to begin to classify and predict how chemicals, like PAHs, produce toxic responses. The developmental toxicity and functional enrichment analyses of the 16 PAHs screened will be addressed in a forthcoming study from the Tanguay laboratory.
Our current study includes several limitations. For example, we used morpholinos, which are a transient method of knocking down gene expression. Another limitation of our study is that we were unable to maintain adequate slincR repression at 72 hpf, which prevents us from fully understanding the functional impact of slincR expression on TCDD-induced toxicity phenotypes that occur later in development. In future studies, we will use a CRISPR/Cas9-generated slincR knockout line to determine the phenotypic impact of slincR expression on later developmental stages. Additionally, our RNA sequencing analysis resulted in a dramatic decrease in the number of significant genes in response to TCDD exposure in the slincR morphants, in comparison with control morphants. One notable limitation of RNA sequencing is the inability to distinguish direct and indirect changes in gene expression; therefore, we are unable to determine if the difference in gene expression is because of the direct regulation of slincR on sox9b and additional target genes or if the difference is because of indirect regulation downstream of slincR target genes. Future studies will identify the genome-wide binding sites of slincR using CHART-sequencing, which has previously shown that lncRNAs can bind hundreds of regions throughout the genome to regulate transcription (West et al. 2014).

Conclusion
We further elucidated the mechanism of TCDD-induced repression of sox9b in the developing zebrafish by showing that slincR is enriched at the 5' UTR of the sox9b locus to repress transcription. We used transcriptome profiling and functional enrichment analyses to demonstrate that ahr2-dependent expression of slincR is involved in processes such as angiogenesis and cartilage and vasculature development. Our data suggest that slincR expression may play a causal role in the TCDD-induced hemorrhaging phenotype, can regulate cartilage development, and is up-regulated by multiple environmentally relevant PAHs. We identified potential mouse and human orthologs of slincR in support of the human health relevance of the zebrafish model. In summary, these data enhance our mechanistic understanding of how AHR activation by environmental pollutants can lead to adverse health effects.