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. 

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.