Chemical modifications of proteins, DNA and RNA moieties play critical roles in regulating gene expression. Emerging evidence suggests the RNA modifications (epitranscriptomics) have substantive roles in basic biological processes. One of the most common modifications in mRNA and noncoding RNAs is N⁶-methyladenosine (m6A). In a subset of mRNAs, m⁶A sites are preferentially enriched near stop codons, in 3′ UTRs, and within exons, suggesting an important role in the regulation of mRNA processing and function including alternative splicing and gene expression. Very little is known about the effect of environmental chemical exposure on m6A modifications. As many of the commonly occurring environmental contaminants alter gene expression profiles and have detrimental effects on physiological processes, it is important to understand the effects of exposure on this important layer of gene regulation. Hence, the objective of this study was to characterize the acute effects of developmental exposure to PCB126, an environmentally relevant dioxin-like PCB, on m6A methylation patterns. We exposed zebrafish embryos to PCB126 for 6 hours starting from 72 hours post-fertilization and profiled m6A RNA using methylated RNA immunoprecipitation followed by sequencing (MeRIP-seq). Our analysis revealed 117 and 217 m6A peaks in the DMSO and PCB126 samples (FDR 5%), respectively. The majority of the peaks were preferentially located around the 3’UTR and stop codons. Statistical analysis revealed 15 m6A marked transcripts to be differentially methylated by PCB126 exposure. These include transcripts that are known to be activated by AHR agonists (e.g., ahrra, tiparp, nfe2l2b) as well as others that are important for normal development (vgf, cebpd, sned1). These results suggest that environmental chemicals such as dioxin-like PCBs could affect developmental gene expression patterns by altering m6A levels. Further studies are necessary to understand the functional consequences of exposure-associated alterations in m6A levels.

Zebrafish embryos were exposed to PCB126 for 6 hours starting from 72 hours post-fertilization and profiled m6A RNA using Methylated RNA immunoprecipitation followed by sequencing (MeRIP-seq). m6A and input RNA libraries were prepared using the SMARTer stranded total RNA-seq kit (Clontech) and sequenced on the Illumina NextSeq500 platform (75bp, paired end reads).

MeRIP sequencing (MeRIP-seq) data analysis involved the following 4 steps.

1. Mapping: Raw reads from individual samples were pre-processed using Trimmomatic to remove adapters and low quality reads. Ribosomal reads were removed from downstream processing by mapping the pre-processed reads to the ribosomal sequences. The zebrafish ribosomal sequences were downloaded from the Silva database (https://www.arb-silva.de/).  Ribosomal-free reads were mapped to the zebrafish genome (GRCz11 version 96) using the STAR aligner version 2.6.0d (Dobin et al. 2013). The resulting unique reads were used for peak calling.

2. Peak calling and visualization: We used Model-based Analysis of ChIP-Seq (MACS) algorithm (version 1.4.2) for identification of the m6A peaks. The immunoprecipitated (m6A) and input (fragmented RNA) samples were used as treatment and control input files, respectively. The effective genome size (--gsize) was set to 117,608,789. The coverage data (wiggle files) was normalized to obtain comparable read density between input and IP samples and visualized using the UCSC genome browser.

3. Motif search and peak annotation: De novo motif finding was done using MEME (Bailey et al. 2009). The FASTA file of the 50 bp region flanking the peak summits was used as input and the top 3 consensus motifs were retrieved. CentriMO was used for visualization of the positional distribution of the best match of the m6A motif in the 300 bp centered around the peak summits.

4. Statistical analysis: Differential methylation analysis was done using DiffBind, a bioconductor package, for differential m6A peak calling (Stark and Brown 2011).