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 AHR-dependent 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 (Garcia et al. 2016; 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 (Garcia et al. 2017). 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 ( and ) to demonstrate that ahr2 is required for the TCDD-induced increase in slincR expression (Garcia et al. 2017; 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 (Garcia et al. 2017). 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 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 control and slincR morphants treated with 0.1% DMSO or 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 . 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 () 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 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 () embryos were exposed to 0, 0.0625, 0.125, 0.25, 0.5, or TCDD (, 95.3% purity; SUPELCO Solutions Within) or vehicle (0.1% DMSO) with gentle rocking for 1 h in glass vials (). The 1-h duration of TCDD exposure was previously shown sufficient to produce the full TCDD-induced toxicity phenotypes (Henry et al. 1997). The concentration was selected to match our previously published TCDD exposure paradigm and results in 99–100% of zebrafish displaying the expected TCDD-induced toxicity phenotype (e.g., cartilage and heart malformations and edema) and an approximate of in sox9b expression (Garcia et al. 2017; 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 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 Petri dishes (at 28°C) in approximately of embryo media until collection at or , depending on the experimental assay. Embryo media consisted of NaCl, KCl, , , , and (Westerfield 2000). The number of embryos raised in a Petri dish ranged from 10–50, depending on the subsequent assay.

Morpholino Injections

Expression levels were reduced using previously published splice-blocking morpholinos (MO) for sox9b (5′-TGCAGTAATTTACCGGAGTGTTCTC-3′), slincR (5′-GACCTAAACTCGACCTTACCAGATC-3′), and a standard negative control (5′-CCTCTTACCTCAGTTACAATTTATA-3′) obtained from Gene Tools, LLC (Garcia et al. 2017; Xiong et al. 2008; Yan et al. 2005). The concentration of slincR morpholino used was demonstrated to produce no visible morphological effects and to knockdown slincR expression at and upon exposure to 0.1% DMSO or TCDD, respectively (Garcia et al. 2017). Approximately, of a [RNA-seq and quantitative real time polymerase chain reaction (qRT-PCR)] or (cartilage morphometrics and hemorrhage screen) solution of slincR MO and control MO were microinjected into the yolk at the one-cell stage. Approximately of a solution of sox9b MO was microinjected into the yolk at the one-cell stage.

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 TCDD, as described above. At , embryos were dechorionated with pronase (). At , embryos were screened for malformations, and the morphologically normal animals were euthanized by overdose of buffered tricaine methanesulfonate (MS-222, ) by prolonged immersion and monitored under a dissecting microscope until their hearts stopped beating (approximately 15 min). Whole 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, of a biotin-labeled capture probe targeting exon 1 of slincR was added to of the sonication solution and incubated overnight at room temperature on an end-over-end 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, 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 zebrafish. The qRT-PCR data were first normalized to 1% input control, such that 6.644 cycles (i.e. dilution factor ) 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 for the two probe sets (). For DNA fold enrichment, we next adjusted relative to the sense probe (), and then fold enrichment was calculated (). The RNA yield was calculated using the following equation (). 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 whole embryos using RNAzol® (Molecular Research Center, Inc.) and a bullet blender with 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 () 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 reactions consisted of 2X SYBR® Green Master Mix (Applied Biosystems), each of forward and reverse primers, and cDNA. Expression values were normalized to and analyzed with the 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 zebrafish with 4 biological replicates per condition (). 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 ) at the shield stage as described in the “Waterborne Exposure” section of the “Methods.” The one-step qRT-PCR reactions were set up consisting of SYBR® Green Master Mix and reverse transcriptase enzyme mix (Power SYBR® Green RNA-to-CT™ 1-Step Kit; Applied Biosystems), each of forward and reverse primers, and 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 and analyzed with the method as described in Livak and Schmittgen (2001). The TCDD concentration served as the calibrator.

Each sample consisted of RNA from 20 pooled 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 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 (, ) 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 .

Zebrafish embryos were exposed to six concentrations of TCDD (0, 0.0625, 0.125, 0.25, 0.5, ) 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 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 as described in Truong et al. (2014).

For each concentration, across two 96-well plates (16 animals per plate per treatment group). The 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 ().

larval photomotor response (LPR) assay.

Zebrafish embryos were exposed to 0, 0.0625, or of TCDD at the shield stage (), as described in the “Waterborne exposure” section of the Methods. At , 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.

For each concentration, 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 , and was compared with the TCDD concentration using a Kolmogorov-Smirnov test ().

PART II: slincR Contributions to TCDD Toxicity

TCDD RNA-seq sample preparation and analysis.

RNA was extracted from control and slincR morphants exposed to 0.1% DMSO or TCDD at with 4 replicates per treatment group, as described above. RNA integrity was confirmed () with an Agilent Bioanalyzer 2,100. Total RNA samples were sent to the Oregon State University Center for Genome Research and Biocomputing Core facilities for library preparation and sequencing. In addition, mRNA was poly A selected, and libraries were prepared with the PrepX™ mRNA and Illumina sequencing workflow (Wafergen Biosystems). Paired-end sequencing () was conducted with an Illumina HiSeq® 3000 sequencer. Processing of next-generation sequencing (NGS) followed a pipeline for gene-level analysis (Anders et al. 2013). Briefly, reads were initially evaluated by FastQC (version 0.11.3) to detect major sequencing problems and quality control trimmed with Skewer (version 0.1.127) to remove ends of reads with low mean Phred quality score using options: –mode pe –end-quality 30 –mean-quality 25 –format auto –threads 9 (Andrews 2015; Jiang et al. 2014). RNA-seq alignment and quantification proceeded with Bowtie2 (version 2.2.3) being used to build TopHat genome index files from the Genome Reference Consortium Zebrafish Build 10 (GRCz10) genome downloaded from Ensembl (release 89), ftp://ftp.ensembl.org/pub/release-89/fasta/danio_rerio/dna/Danio_rerio.GRCz10.dna.toplevel.fa.gz (Kim et al. 2013; Langmead et al. 2009; Langmead and Salzberg 2012). Trimmed reads were aligned with TopHat (version 2.1.1) using the following options optimized for trimmed, paired-end, stranded data: –num-threads 9 –library-type fr-secondstrand –no-discordant –no-mixed –mate-inner-dist 20 –mate-std-dev 180 –min-anchor-length 5 –min-intron-length 10 –max-intron-length 400,000. Samtools (version 1.4) was used to assess the percent of paired aligned reads and to sort aligned reads (binary alignment files, BAM) by name using the “flagstat” and “sort -n -@9” commands, respectively (Li et al. 2009). Gene counts were estimated using the htseq-count script from HTSeq (version 0.6.0) with the GRCz10 Ensembl (release 89) GTF annotation, ftp://ftp.ensembl.org/pub/release-89/gtf/danio_rerio/Danio_rerio.GRCz10.89.gtf.gz, and options: –format=bam –idattr=gene_id –stranded=yes –mode=intersection-nonempty –order=name (Anders et al. 2015).

Each biological sample consisted of RNA from 20 pooled zebrafish with 4 biological replicates per condition (). Differential expression analysis followed the maintained Bioconductor workflow developed by Y Chen et al. (2016), 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; (Huber et al. 2015). The Bioconductor package edgeR (version 3.22.3) was used to normalize counts and identify differentially expressed genes (Lun et al. 2016; McCarthy et al. 2012; Robinson and Smyth 2007, 2008; Robinson et al. 2010; 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 Chen et al. 2016; 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 “” option was used to protect the empirical Bayes estimates against the possibility of outlier genes with wide-ranging individual dispersions. Genes with a Benjamini-Hochberg (BH) adjusted 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 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 TCDD, as described above. At , larvae were euthanized by overdose of buffered tricaine methanesulfonate (MS-222, ) 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% . The cartilage was stained with 0.4% Alcian Blue 8GX (Sigma-Aldrich) in 70% ethanol and 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 () 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).

Biological replicates consisted of 9–10 individual zebrafish measurements per treatment group. Statistical significance ( or 10, ) was determined in RStudio (version 1.1.453) with the R (version 3.5.1; R Development Core Team) packages multcomp (version 1.4-8) and sandwich (version 2.4-0) using a modified one-way ANOVA and Tukey multiple comparison of means that is robust in regard to equal variance, sample sizes, and distribution, and data were graphed using ggplot2 (version 3.0.0), and a custom R script (Cartilage_morphometrics.R; Supplemental Material, S5. Cartilage morphometrics R script; Herberich et al. 2010; Hothorn et al. 2008; Wickham 2016; Zeileis 2004, 2006).

Hemorrhage screen.

Zebrafish embryos were microinjected with control and slincR MOs and exposed to 0.1% DMSO or TCDD, as described above. At , the embryos were evaluated under a dissecting microscope for the presence or absence of hemorrhaging.

The DMSO and TCDD samples consisted of 3 or 5 biological replicates, respectively. Each biological replicate contained 10–12 individual 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 ( or 5, ) 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 zebrafish using a custom automated dechorionator. At , 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 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 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 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 , ). Of the 16 PAHs, 9 did not result in morphological malformations and were subsequently tested at only . The remaining 7 were tested using a definitive concentration range (11 concentrations, ) 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 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 log-logistic 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 (). The computed 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).

Table 1.

RNA sequencing gene expression results from zebrafish developmentally exposed to 16 PAHs from the shield stage to .

Treatment	cyp1a		slincR		sox9b
	log_2FC	FDR	log_2FC	FDR	log_2FC	FDR
Retene (12.2μM)	7.692	6.94×10^−16	3.006	2.50×10^−10	0.071	0.848
Benzo[j]fluoranthene (14.9μM)	7.761	4.96×10^−16	3.922	9.65×10^−14	0.095	0.678
Benzo[k]fluoranthenea (1.9μM)	7.913	2.91×10^−16	3.975	6.07×10^−14	0.168	0.792
Dibenzo[a,h]pyrenea (5μM)	6.373	1.07×10^−13	2.936	1.09×10^−9	0.104	0.610
Dibenzo[a,i]pyrene (5μM)	6.478	7.58×10^−14	2.785	5.17×10^−9	−0.191	0.811
Benzo[b]fluoranthenea (50μM)	6.019	1.83×10^−13	2.139	4.25×10^−6	0.171	0.615
Fluoranthenea (50μM)	3.142	1.40×10^−6	0.343	0.720	0.001	0.998
Phenanthrenea (50μM)	2.093	0.0026	0.009	1.00	0.011	1.00
Acenapthenea (50μM)	0.844	0.139	0.459	0.407	0.113	0.499
4H-cyclopenta[def]phenanthren-4-one (16.2μM)	1.618	0.176	0.019	0.990	0.347	0.0664
Carbazole (50μM)	2.126	0.0549	0.491	0.716	0.211	0.384
3-nitrofluoranthene (1.9μM)	1.778	0.712	0.319	0.954	0.186	0.859
1,5-dimethylnaphthalene (50μM)	0.806	1.00	0.401	1.00	0.003	1.00
9-methylanthracene (50μM)	2.647	0.0091	0.711	0.411	0.140	0.511
2-methylnaphthalene (50μM)	0.587	1.00	0.139	1.00	0.038	1.00
Anthracenea (50μM)	2.467	0.366	1.083	0.415	0.137	0.632

Note: These data are part of an unpublished RNA sequencing dataset that was generated by Dr. M. Geier from the R. Tanguay laboratory. Concentrations were based on the and were experientially determined by the individual compound response. For the PAHs that did not elicit a morphological effect, the maximum soluble concentration tested ( for most) was used. The Bioconductor package edgeR (statistical methodology based on the negative binomial distribution) was used to calculate the and identify the significant differentially expressed genes (, BH-adjusted ), when compared to the vehicle control (0.1% DMSO). ; ; ; ; ; .

Denotes a chemical that is one of the U.S. Environmental Protection Agency’s 16 priority PAHs.

Total RNA was isolated and quantified from pooled groups of eight 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 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.

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 zebrafish with 4 biological replicates per condition (, BH-adjusted ).

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 TCDD on E13.5 following previously established protocols (Lin et al. 2003). At E16.75, dams were euthanized via CO2 asphyxiation, and the urogenital sinuses (UGS) were dissected from fetuses as previously described (Branam et al. 2013). Each UGS was immediately placed into a microfuge tube containing of 1% trypsin (Difco, 215240) in PBS and incubated on ice for 30 min. Collagenase (Sigma C9891) was added to a final concentration of , 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-mixed –num-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.

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 (, BH-adjusted ).

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 (Chen et al. 2016). 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 (BioProject: PRJNA66167; Gene ID: 100499467) and SOX9 (BioProject PRJEB4337; Gene ID: 6662). The data were graphed using the R package ggplot2 (version 3.0.0; Fagerberg et al. 2014; Wickham 2016; Yue et al. 2014).

Results

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′ 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 CHART-isolated RNA; however, the RNA samples did not result in detectable amplification (no measurable enrichment).

Figure 1.

Use of qRT-PCR CHART to determine slincR enrichment at the sox9b promoter in whole embryos treated with 0.1% DMSO or TCDD. (A) Enrichment of slincR RNA by CHART (slincR probe) and a nonspecific primer set (). Each condition had 3 biological replicates, where 1 replicate consisted of approximately 500 zebrafish. The qRT-PCR data were first normalized to 1% input control, such that 6.644 cycles (i.e., dilution factor ) 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 for the two probe sets (). The RNA yield was calculated using the following equation (). 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 (, , and ) downstream of the sox9b transcription start site and 5' untranslated region. Each condition had 3 biological replicates (), where 1 replicate consisted of approximately 500 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 (), and then fold enrichment was calculated (). 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 whole embryo control and sox9b morphants. Expression values were analyzed with the method and normalized to , 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 (), and the data were analyzed using a one-way ANOVA with a Tukey post hoc test ( in comparison with ). Error bars (A, C, and D) indicate standard error of the mean. Note: ; ; ; ; .

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 . 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 of TCDD and measured the relative expression levels of sox9b and slincR in 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 TCDD concentration and increased with increasing concentrations of TCDD, in comparison with TCDD (i.e., 0.1% DMSO; Figure 2A). The expression levels of sox9b were significantly decreased in the 0.5 and TCDD concentrations, in comparison with TCDD (Figure 2A).

Figure 2.

Analysis of the concentration–response effects of developmental exposure (1-h exposure at ) to TCDD on gene expression at and morphological malformations at . (A) qRT-PCR relative expression of cyp1a, slincR, and sox9b transcripts in wild-type embryos exposed to 0, 0.625, 0.125, 0.25, 0.5, or TCDD. For all assays, 0.1% DMSO served as the vehicle control and is listed as TCDD. Expression values were analyzed with the method and normalized to using the TCDD concentration as the calibrator. Each experimental unit represents a pool of 20 embryos, and each treatment group included four biological replicates (). The data were and analyzed using a one-way ANOVA with a Dunnett post hoc test ( in comparison with ). All error bars indicate standard error of the mean. (B) Evaluation of 17 physical malformations at on wild-type zebrafish exposed to six concentrations of TCDD () 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 (, ). (C) Representative lateral images of zebrafish for each concentration of TCDD tested. The bar in top left corner indicates . (D) Larval photomotor response (LPR) in wild-type embryos developmentally exposed to 0, 0.0625, or 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 TCDD using a Kolmogorov-Smirnov test (, ). Note: ; ; ; ; ; ; ; ; ; at ; at ; ; ; ; .

In addition to gene expression changes, we also examined the embryos for developmental toxicity when exposed to the six concentrations of TCDD () using a 17 end point morphology screen at , as described in Truong et al. (2014). Developmental exposures of 0.25, 0.5, or 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 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 TCDD (Figure 2D). Larval zebrafish exhibiting any visible malformation were necessarily excluded, which prevented the analysis of the 0.25, 0.5, or TCDD concentrations in the LPR assay.

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 slincR and control morphants exposed to 0.1% DMSO or TCDD. In control morphants exposed to TCDD, we identified 132 genes that were differentially expressed when compared with vehicle control (, BH-adjusted , Figure 3A and Excel Table S1). In slincR morphants exposed to TCDD, we identified 80 differentially expressed genes when compared with vehicle control (, BH-adjusted , 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).

Figure 3.

Venn diagram of all differentially expressed genes, and heatmaps of significantly differentially expressed genes with the top 30 largest between slincR and control morphants treated developmentally (1 h exposure at ) with 0.1% DMSO or TCDD. The Bioconductor package edgeR (statistical methodology based on the negative binomial distribution) was used to identify the significant differentially expressed genes in TCDD-treated morphants (, BH-adjusted ), in comparison with the vehicle control (0.1% DMSO). (A) Venn diagram of up- and down regulated genes from control and slincR morphants. Hierarchical clustering is derived from TMM-normalized, regular log-transformed gene counts scaled by z-score as shown by heatmaps of (B) control morphant TCDD (CT) in comparison with control morphant DMSO (CD) and (C) slincR morphant TCDD (ST) in comparison with slincR morphant DMSO (SD). slincR is represented by the Ensembl gene symbol FP017215.1 in the heatmaps. Each condition included four biological replicates (i.e., CD1, CD2, CD3, and CD4), with each replicate representing a pool of 20 embryos. Note: ; ; .

To visualize the expression profile, we performed bidirectional hierarchical clustering of the differentially expressed genes with the top 30 largest 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 ) 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).

Table 2.

Significantly enriched MetaCore process networks from Con-MO zebrafish developmentally exposed to TCDD for 1 h at the shield stage, compared to the vehicle control (0.1% DMSO).

Network name	FDR	Associated human gene
Development: Cartilage development	8.053×10^−3	BMP4, BMP6, COL12A1, SOX9, CTGF
Cell adhesion: Cell-matrix interactions	8.053×10^−3	COL12A1, THBS1, MMP9, VCAN, ITGA11, TNC, COL10A1
Immune response: Th17-derived cytokines	8.053×10^−3	CXCL12, IL17RC, MMP9, JUN
Development: Ossification and bone remodeling	8.053×10^−3	BMP4, BMP6, OSTN, FOXO1, COL10A1, POSTN
Signal Transduction: BMP and GDF signaling	8.053×10^−3	BMP4, BMP6, CDKN1B, GADD45B, SOX9
Development: EMT: Regulation of epithelial-to-mesenchymal transition	8.053×10^−3	BMP4, HMGA2, MMP9, JUN, CDKN1B, SOX9, CTGF
Proteolysis: ECM remodeling	8.053×10^−3	MMP9, TNC, COL10A1, CTGF
Signal Transduction: ESR1-nuclear pathway	9.733×10^−3	CYP1A1, CXCL12, BMP6, ADM, JUN, FOXO1, CYP1B1
Reproduction: FSH-beta signaling pathway	2.218×10^−2	BMP4, JUN, FOXO1, CDKN1B, CTGF
Development: Blood vessel morphogenesis	2.424×10^−2	EPAS1, CXCL12, BMP4, FOXO1, CTGF, ANGPT4

Note: The Bioconductor package edgeR (statistical methodology based on the negative binomial distribution) was used to calculate the significant differentially expressed genes (, BH-adjusted ), 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 . ; ; ; ; ; ; ; ; ; ; ; .

Table 3.

Significantly enriched MetaCore process networks from slincR-MO zebrafish exposed to TCDD for 1 hour at the shield stage, in comparison with the vehicle control (0.1% DMSO).

Network name	FDR	Associated human gene
Immune response: Th17-derived cytokines	2.946×10^−2	JUN, CXCL12, MMP9
Development: EMT Regulation of epithelial-to-mesenchymal transition	2.946×10^−2	JUN, CDKN1B, TWIST1, CTGF, MMP9
Signal transduction: ESR1-nuclear pathway	2.946×10^−2	CYP1A1, JUN, NFATC4, CXCL12, CYP1B1
Reproduction: FSH-beta signaling pathway	2.946×10^−2	JUN, CDKN1B, FST, CTGF

Note: The Bioconductor package edgeR (statistical methodology based on the negative binomial distribution) was used to calculate the significant differentially expressed genes (, BH-adjusted ), 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 . ; ; ; ; ; ; ; .

Table 4.

Significantly enriched unique MetaCore GO processes related to skeletal and cartilage development from Con-MO zebrafish exposed to TCDD for 1 h at the shield stage, in comparison with the vehicle control (0.1% DMSO).

Biological process	FDR	Associated human gene
Ossification	7.152×10^−5	CYP24A1, BMP4, BMP6, MMP9, OSTN, VCAN, ITGA11, TNC, SOX9, COL10A1, CTGF
Cartilage development	2.601×10^−4	BMP4, BMP6, HMGA2, THBS1, SFRP2, SOX9, COL10A1, CTGF
Endochondral bone morphogenesis	1.968×10^−3	BMP4, BMP6, THBS1, SOX9, COL10A1
Osteoblast differentiation	4.093×10^−3	CYP24A1, BMP4, BMP6, VCAN, ITGA11, TNC
Bone development	5.122×10^−3	BMP4, BMP6, THBS1, PAPPA2, SFRP2, SOX9, COL10A1
Regulation of chondrocyte differentiation	5.710×10^−3	BMP4, BMP6, SOX9, CTGF
Chondrocyte differentiation	1.142×10^−2	BMP4, HMGA2, SFRP2, SOX9
Regulation of cartilage development	1.177×10^−2	BMP4, BMP6, SOX9, CTGF
Positive regulation of osteoblast differentiation	1.271×10^−2	IL6R, BMP4, BMP6, SFRP2
Regulation of osteoblast differentiation	1.324×10^−2	IL6R, BMP4, BMP6, OSTN, SFRP2
Regulation of bone mineralization	1.356×10^−2	BMP4, BMP6, OSTN, SOX9
Odontogenesis of dentin-containing tooth	1.356×10^−2	BMP4, ADM, FOXO1, TNC
Positive regulation of cartilage development	1.567×10^−2	BMP4, BMP6, SOX9
Regulation of biomineral tissue development	1.578×10^−2	BMP4, BMP6, OSTN, SOX9
Regulation of ossification	1.793×10^−2	IL6R, BMP4, BMP6, OSTN, SFRP2, SOX9
Chondrocyte proliferation	1.793×10^−2	HMGA2, CTGF
Endochondral ossification	2.003×10^−2	BMP4, BMP6, COL10A1
Replacement ossification	2.003×10^−2	BMP4, BMP6, COL10A1
Positive regulation of ossification	2.891×10^−2	IL6R, BMP4, BMP6, SFRP2
Odontogenesis	2.954×10^−2	BMP4, ADM, FOXO1, TNC
Negative regulation of chondrocyte differentiation	4.076×10^−2	BMP4, SOX9

Note: The Bioconductor package edgeR (statistical methodology based on the negative binomial distribution) was used to calculate the significant differentially expressed genes (, BH-adjusted ), 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 . ; ; ; ; ; ; .

Table 5.

Significantly-enriched unique MetaCore GO processes related to angiogenesis and vasculature development from slincR-MO zebrafish exposed to TCDD for 1 h at the shield stage, in comparison with the vehicle control (0.1% DMSO).

Biological process	FDR	Associated human gene
Regulation of blood vessel endothelial cell proliferation involved in sprouting angiogenesis	9.926×10^−3	HMOX1, THBS1
Negative regulation of vascular smooth muscle cell proliferation	1.297×10^−2	CDKN1B, HMOX1
Regulation of blood vessel endothelial cell migration	2.806×10^−2	HMOX1, NFE2L2, THBS1
Positive regulation of vascular smooth muscle cell proliferation	2.941×10^−2	JUN, MMP9
Positive regulation of blood coagulation	3.062×10^−2	NFE2L2, THBS1
Regulation of cell migration involved in sprouting angiogenesis	4.019×10^−2	HMOX1, THBS1
Branching involved in blood vessel morphogenesis	4.456×10^−2	NFATC4, CXCL12

Note: The Bioconductor package edgeR (statistical methodology based on the negative binomial distribution) was used to calculate the significant differentially expressed genes (, BH-adjusted ), 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 compared to the background gene list. Enriched GO processes were considered significant with an . ; ; ; ; ; ; .

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 TCDD for 1 h at the shield stage () and evaluated the cartilage of zebrafish larvae. Both the control and slincR morphants exposed to TCDD had a 100% incidence of cartilage malformations at . The cartilage of DMSO-exposed 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 ( or 10, ; 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).

Figure 4.

Effect of TCDD on the developing jaw after treatment of control and slincR morphants. Control (Con-MO) and slincR morphants (slincR-MO) were developmentally exposed to 0.1% DMSO or TCDD, and the cartilage was stained and measured at . (A) A morphometric system was used to measure the position and length of landmark structures in the developing jaw. The position of jaw structures representing junctions between Meckel’s and palatoquadrate cartilages (point A and B) and representing junctions between the hyosymplectic and ceratohyal cartilages (point C and D) was measured relative to a reference point as shown. Statistical significance was determined by a modified two-way ANOVA with a Tukey post hoc test. Morphometric values represent . ( individual zebrafish; , , ). The asterisk (*) indicates statistical significance between slincR-MO and Con-MO samples exposed to TCDD. (B) Representative images of the cartilage structures in control and slincR morphants treated with DMSO or TCDD. The arrow in the TCDD Trace panel points to the junction between the hyosymplectic and ceratohyal cartilages (point C and D). The bar in the top right corner indicates . Note: ; ; .

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 TCDD for 1 h at the shield stage (). At , 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 ( or 5, ; 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).

Figure 5.

TCDD-induced hemorrhaging in control and slincR morphants. Control (Con-MO) and slincR morphants (slincR-MO) were developmentally exposed to 0.1% DMSO or TCDD, and 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 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 ( or 5, ). (B) Representative images of wild-type zebrafish exposed to DMSO (normal phenotype) or TCDD (hemorrhaging phenotype). The bar in the top right corner indicates and arrows point to sites of hemorrhaging. Note: ; ; .

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 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.

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 of TCDD, 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; Chen et al. 2016). A sequence comparison of the zebrafish slincR sequence () with the mouse () and human () 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.

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 ( and ). 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: ; .

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 (Noyes et al. 2016). We previously used the zebrafish model to identify an ahr2-dependent long noncoding RNA, slincR, whose expression increased upon exposure to multiple Ahr2 ligands (Garcia et al. 2017). 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 (Garcia et al. 2017).

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 (Garcia et al. 2017). 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 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 ) were able to both significantly increase slincR and significantly decrease sox9b expression levels. In response to developmental exposure to TCDD, sox9b expression is reduced in the eye, brain, otic vesicle, jaw region, and snout region (Garcia et al. 2017). For the TCDD concentration, the discrepancy 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.

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 exposed to 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 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 in head mesenchyme provided in vivo evidence that suggests 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 mutant displayed a significant increase in expression of Sox9 and its upstream regulator Runx2, which argues that is not directly regulating Sox9. In zebrafish, several publications have supported the hypothesis that xenobiotic activation of the AHR leads to a disruption in 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 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 signaling. In addition, exposure to a pharmacological agent (6-bromoindirubin-3′-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 ahr2-dependent expression of slincR plays a causal role in TCDD-induced 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 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 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 (Garcia et al. 2017). 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 , 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.

Click here for additional data file.

Click here for additional data file.

Click here for additional data file.