Whole genome bisulfite sequencing (WGBS) analysis

WGBS libraries were prepared using the TruSeq DNA Methylation Library Preparation Kit (Illumina). Briefly, genomic DNA was isolated using the DNeasy Blood & Tissue Kit (Qiagen) according to the manufacturer’s guide. Then, approximately 500 ng DNA were bisulfite treated using EZ DNA Methylation Kit (Zymo Research). Bisulfite-converted DNA was end-repaired, dA-tailed, and ligated with adapters following instructions of the TruSeq DNA Methylation Library Preparation Kit. Finally, high quality indexed libraries were then pooled and sequenced on Illumina NovaSeq platform with 150-bp paired-end reads.

WGBS data analysis was followed our established analysis pipelines.^57,58 Briefly, WGBS raw data were removed first 12-bp at the 5’ end of both pairs, and reads with adapters and low-quality bases by using TrimGalore-0.4.3. The trimmed sequences were mapped to the bovine genome (ARS-UCD1.2) using Bismark. Uniquely mapped reads were then removed PCR duplicated reads and non-converted reads using deduplicate_bismark and filter_non_conversion. For avoiding the sequencing bias, only reads with 10x coverage was used in the downstream analysis. Methylation of each CpG site was calculated and methylation DNA methylation of each sample was calculated by averaging the consecutive genomic window of 300-bp tiles’ methylation. Differentially methylated regions (DMRs) were defined as common 300-bp tiles between two compared groups, which methylation levels ≥ 75% in one group, while ≤25% in another, and were significantly different by Fisher’s exact test (P-value ≤0.05, FDR ≤0.05). Hyper- and hypo-methylated tiles were those with DNA methylation levels ≥75% and ≤25%, respectively. The gene ontology and pathway analysis were performed by means of the David tool.⁵⁶