meDIP-seq methods

Bioinformatic analysis of meDIP-seq samples was performed by aligning paired-end sequenced FASTQs to the porcine reference genome (susScr11.1; using bowtie2 (v 2.3.3.1) [31]. Alignment files were filtered to retain properly mapped pairs in which at least one-end maps uniquely to the genome. PCR duplicates were removed using PICARD (v1.67) MarkDuplicates. The meDIP peaks were called using MACS [32] to identify punctate binding events with both high sensitivity and specificity [33].

Differential Peak Analysis was performed using the DiffBind (v2.14.0) [34] application package. This package executes differential binding analysis of ChIP-Seq, and similar epigenomics profiling techniques, peak data by comparing biological replicates from different experimental conditions and determines sites of differential 5hmC coverage. DiffBind analysis was performed on peaks that were called to be enriched by MACS, as well as the genomic coordinates for exons and promoters in Ensembl (release 98) [35]). Differential peaks and genomic coverage bins were annotated and assigned to the corresponding genes by the HOMER (v4.10) [36] peak annotation tool.

Manual analysis of 5hmC coverage was performed using per-base coverage of regions of interest, calculated by bedtools (v2.20.0) genomeCoverageBed. Sequence read values for the overall exonic 5hmC coverage per gene were also calculated in parallel using htseq-count (v0.9.1) [37]. These faux-RNA counts 5hmC coverage were then processed by edgeR (v3.28.1) [38] to obtain differences in read frequencies for genomic 5hmC coverage at exons analogous to differential expression analysis of RNA reads.