DOI: 10.5061/dryad.7pvmcvdw1
title: Exposure to 3, 3’, 4, 4’, 5-pentachlorobiphenyl (PCB 126) causes widespread DNA hypomethylation in adult zebrafish testis

description: Exposure to environmental toxicants during preconception have been shown
 to affect offspring health and epigenetic mechanisms such as DNA
 methylation are hypothesized to be involved in adverse outcomes. However,
 studies addressing the effects of exposure to environmental toxicants
 during preconception on epigenetic changes in gametes are limited. The
 objective of this study is to determine the effect of preconceptional
 exposure to a dioxin-like PCB (PCB126) on DNA methylation and gene
 expression in testis. Adult zebrafish were exposed to 3 and 10 nM PCB126
 for 24 hours and testis tissue was sampled at 7 days post-exposure for
 histology, DNA methylation and gene expression profiling. Reduced
 Representation Bisulfite Sequencing revealed 37 and 92 differentially
 methylated regions (DMRs) in response to 3 and 10nM PCB126 exposures,
 respectively. Among them 19 DMRs were found to be common between both
 PCB126 treatment groups. Gene ontology (GO) analysis of DMRs revealed that
 enrichment of terms such as RNA processing, iron-sulfur cluster assembly
 and gluconeogenesis. Gene expression profiling showed differential
 expression of 40 and 1621 genes in response to 3 and 10nM PCB126
 exposures, respectively. GO analysis revealed differential expression of
 genes related to xenobiotic metabolism, oxidative stress and immune
 function. There is no overlap in the GO terms or individual genes between
 DNA methylation and RNAseq results, but functionally many of the altered
 pathways have been shown to cause spermatogenic defects. Our results
 indicate that exposure to dioxin-like PCBs during preconception could
 affect testicular function by altering DNA methylation patterns, with
 significant implications for reproductive health. 
description: This data is part of an experimental study investigating the
 effect of exposure to PCB126, a dioxin-like PCB congener on DNA
 methylation and gene expression profiles in zebrafish testis. Adult male
 zebrafish (6 months old) were exposed to either 3 or 10 nM PCB126 or
 solvent carrier (0.01% DMSO) for 24 h. At the end of the exposure, fish
 were transferred to clean water with constant aeration and heating and
 maintained for 7 days. Fish were not fed during the 24 h exposure period.
 At 7 days post-exposure, fish were euthanized with MS-222 (150 mg/l)
 buffered with sodium bicarbonate (pH 7.2) prior to tissue sampling. We
 chose this experimental design in order to capture both primary and
 secondary changes in DNA methylation and gene expression.
 DNA methylation profiling was done by Reduced Representation
 Bisulfite Sequencing (RRBS). RRBS library preparation and sequencing was
 conducted by NXT-Dx, a Diagenode company (Ghent, Belgium).  Raw reads were
 pre-processed using TrimGalore (v0.4.3) and cutadapt (v1.12). FASTQC
 (v.0.10.1) was used to determine the quality of the sequencing reads.
 Pre-processed reads were aligned to the zebrafish genome (GRCz10) using
 the Bisulfite Analysis Toolkit (BAT; (Kretzmer et al. 2017)). BAT is an
 integrated toolkit that includes aligning the reads to the genome, calling
 differentially methylated regions (DMRs) using metilene (Juhling et al.
 2016), annotation of DMRs, statistical analysis and correlating DMRs with
 gene expression. This is in contrast to prevalent methods where separate
 pipelines are used for mapping, DMR calling and statistical analysis. A
 comparison of different DNA methylation profiling pipelines to identify
 DMRs using RRBS data demonstrated that in most cases metilene has the
 highest precision (Liu et al. 2020). After the read
 mapping the sequencing runs, one control (D1) and one 3nM PCB126 group
 (P3-12) were considered as outliers and were excluded from the subsequent
 analysis. Only cytosine positions in a CpG context with at least 10 and at
 most 100 overlapping reads were considered. The DMRs were required to
 contain at least 10 of such cytosines. The minimum difference of mean
 methylation rates per group was at least 0.1 and only DMRs with a q-value
 of 0.05 were considered significant. We classified the DMRs into various
 genic (promoters, introns and exons) and intergenic regions, CpG island
 and shores and repetitive elements using UCSC table browser.
 We used Genomic Regions Enrichment of Annotations Tool (GREAT) to
 associate DMRs with genes (McLean et al. 2010). GREAT predicts gene
 functions of cis regulatory elements by assigning each gene a regulatory
 domain. To use GREAT, we converted the genomic coordinates of DMRs from
 GRCz10 version to Zv9 version of the genome using the UCSC genome browser
 liftOver utility (https://genome.ucsc.edu/cgi-bin/hgLiftOver). We used default parameters with a basal domain that extends 5 kb upstream and 1 kb downstream of the TSS and conducted gene ontology (GO) (biological process and molecular function) enrichment analysis. Stranded RNAseq library preparation using the Illumina TruSeq total RNA library prep kit and 50 bp single-ends sequencing on the HiSeq2000 platform were performed at the Tufts University Core Facility. Raw data files were assessed for quality using FastQC  (Andrews 2010) and pre-processed as described previously. Data analysis was done as described previously (Aluru et al., 2018) by mapping the pre-processed reads to the Ensembl version 90 (GRCz10) of the zebrafish genome. Mapped reads were counted using HTSeq-count (Anders et al. 2015). Statistical analysis was conducted using DESeq2, a Bioconductor package (Love et al. 2014). Only genes with false discovery rate (FDR) of less than 5% were considered to be differentially expressed. Gene Ontology (GO) term enrichment analysis was performed using gProfiler package g:GOSt (Reimand et al. 2016). The up- and down-regulated differentially expressed genes (DEGs) from the high PCB126 group were analyzed separately. Ensembl Gene IDs of DEGs were used as input. The resulting output of significantly enriched GO Biological Process (BP) and Molecular Function (MF) terms were used as input in REVIGO (Supek et al. 2011) to remove redundant GO terms. The resulting GO terms were used to draw visualizations using CirGO (Kuznetsova et al. 2019).  CirGO allows hierarchical visualization, where the inner circle represents the parent terms and the outer circle shows descendant child terms. Child terms are sorted based on statistical significance and highlighted with a color gradient. Raw and processed RRBS and RNAseq data have been deposited into Gene Expression Omnibus (Accession number GSE190741). In addition, the RRBS and RNAseq data can be visualized through the UCSC genome browser track hub: http://www.bioinf.uni-leipzig.de/~jane/Neel/testisRRBSHub/hub.txt. The instructions for visualizing the data are provided in the supplementary information.