A total of 5 pure blood Angus bulls were selected for semen collection. Semen samples were collected twice from each bull during early March 2019 designated as Winter samples and early June 2019 designated as Late Spring samples (Supplementary Table S3). Semen was collected from 5 pure blood Angus bulls by ABS Global, DeForest, WI (Supplementary Table S3). Because bovine sperm formation requires 61 days (Staub and Johnson, 2018), bulls were exposed to the coldest temperature (January and February) during spermatogenesis for semen collected in March (Winter group), and warm temperature (April and May) for semen collected in the early June (Late Spring group). During the winter of 2019, the average temperature during January to March was 0.5°C (max 6 to min -5°C) and during late spring the average daily temperature during the April and May were 8.01°C (max. 21 to min. 10°C) and 13.2°C (max. 28 to min. 18°C) respectively. All the semen straws were transported to the lab in a liquid nitrogen shipper and stored at −80°C temperature.

Semen samples from the same 3 randomly selected bulls out of 5 bulls collected during the winter and late spring respectively were used for isolation of genomic DNA as previously described (Perrier et al., 2018). One straw from each bull was used for DNA extraction (about 20 million spermatozoa). After thawing a straw at 37°C temperature, semen extender was removed by Phosphate Buffer Saline (PBS) wash and incubated with lysis buffer (10 mM Tris-HCl pH 7.5, 25 mM EDTA, 50 mM dithiothreitol, 75 mM NaCl, 1% SDS and 0.5 μg glycogen)) and proteinase K (0.2 mg/ml) at 55°C overnight. The phenol: chloroform (1:1) was used to extract the DNA after incubation with RNAse A for 1 h at 37°C and washed with ethanol. Extracted DNA quality was checked by using NanoDrop 2000 spectrophotometer (Nanodrop Technologies, USA).

Bisulfite conversion and WGBS library were prepared by following a previously described method by Novogen (Zhou et al., 2018). High quality genomic DNA (2 µg) from each sample was used for library preparation and bisulfite conversion. Briefly, the genomic DNA spiked with lambda DNA were fragmented to 200–400 bp with a Covaris s220 sonicator (Covaris, Inc, USA) followed by terminal repair and adenylation. Lambda DNA was used as unmethylated control for bisulfite conversion rate calculation. Cytosine methylated barcodes were then ligated to the sonicated DNA and treated twice with an EZ DNA Methylation-Gold™ Kit (Zymo Research) for bisulfite conversion. Then, bisulfite converted libraries were sequenced on an Illumina HiSeq 2000/2,500 platform at the Novogene Bioinformatics Institute (Beijing, China).

FastQC was used to generate quality score of the sequences. Sequences with a PHRED score <20 and all adapter sequences were trimmed with Trim Galore program v 5.0 (Martin, 2011). The clean data from each sample were merged to align with the reference genome. The Cattle reference genome ARS-UCD 1.2 (GCA 002263795.2 https://ftp.ensembl.org/pub/release-108/gtf/bos_taurus/) incorporated with Y chromosome (GenBank: CM011803.1) were used for paired end alignment by indexing with bowtie2 (Version 2.4.5) under bismark (0.22.3, released: 19-11–2019) (Krueger and Andrews, 2011; Langmead and Salzberg, 2012).

Bismark deduplication script “deduplicate_bismark” was used to remove reads in the same location from alignment output and the methylation information was extracted in CpG/CHG/CHH context by using “bismark_methylation_extractor” tool. To avoid the effect of bias towards non-methylation in the end of read due to end repairing, first 6 bp were ignored (Krueger and Andrews, 2011).

Two popular Bioconductor packages Methylkit and DSS of R platform were used to perform the downstream analysis (Akalin et al., 2012; Feng and Wu, 2019). Sequence bases contained less than 10 uniquely mapped reads were filtered and removed across all the samples by Methylkit. Sequence bases containing at least 10 uniquely mapped read were termed as CpG10s and used for further analysis. Methylation percentage in each CpG site was calculated using the DMLtest function of DSS where a single CpG site containing greater than 10% difference and a p-value < 0.05 between the Late Spring and Winter groups were termed as a DMC. A DMC containing a higher mean methylation percentage during Winter in comparison with Late Spring group was termed as Hypermethylated DMC and similarly lower mean methylation percentage in Winter was term as hypomethylated DMC.

DMRs are genomic regions containing several adjacent DMCs and possessing a higher or lower methylation percentage in the experimental group. To form a DMRs, DSS first identifies statistically significant DMCs and merges several adjacent DMCs into a region (Feng and Wu, 2019). We calculated the mean methylation of CpG10s across all samples by using “dmlTest” function of DSS and generated the DMRs using default “callDMR” function with a p-value < 0.05. We also enabled the smoothing option of dmlTest function for besting estimation of mean methylation.

A popular R package Genomation operated on another R package GRanges was used for annotation of DMRs (Lawrence et al., 2013; Akalin et al., 2015). DMRs overlapped with different regions of the genes (promoters, introns, and exons) and distance from TSS were calculated by using “readTranscriptFeatures” and “getAssociationwithTSS” function of Genomation package. Promoters were defined as −2000 to −100 bp relative to the TSS sites. The association of DMRs with CGI was calculated by downloading CG data from California Santa Cruz Table Browser (Karolchik et al., 2004). During annotation with CGI, a CGI shore was defined as 4 kbp distance on the either side of each CGI. The list of annotated DMCs and DMRs is available in the Additional File 2: Supplementary Table S1 and Additional File 3: Supplementary Table S2.

All unique genes containing at least one DMR were further used for Gene function enrichment analysis by Database for Annotation, Visualization and Integrated Discovery (DAVID) (Huang da et al., 2009). Clusters of terms generated by DAVID showing EASE enrichment scores above 0.5 were considered significant. Along with this, Gene Ontology Enrichment analysis was performed by Gene Ontology Consortium’s online tool (http://www.geneontology.org/) to annotate the DMGs in Biological processes, Molecular functions and Cellular component categories (Mi et al., 2019). The identified unique genes were then classified and grouped into different pathways by using Reactome pathway and p-value <0.05 was considered as significant for each group (Fabregat et al., 2018).

To verify the DMRs, we performed MS-PCR of 6 genes Smg9, Mest, Peg10, Lmbr1, C21H15orf40 and Prmt6. Genomic DNA was isolated from sperm samples with Proteinase K digestion. The quality of extracted DNA was checked using NanoDrop 2000 spectrophotometer (Nanodrop Technologies, USA) and samples with a A260/280 score above 1.8 were used for further analysis. Extracted DNA samples were used for Bisulfite conversion by EZ DNA Methylation-Direct™ Kit (Catalog No. D5020, Zymo Research, Irvine, CA, USA). Briefly, around 2 µg of the extracted DNA was mixed with 130 µL bisulfite conversion (BS) reagents and placed in a thermal cycler for bisulfite conversion with the following condition: 98°C for 8 min, 64°C for 3.5 h and stored 4°C until further use. Through BS conversion, unmethylated cytosines (C) were converted to uracil (U) while methylated cytosines remained unchanged. BS converted DNA samples were amplified by PCR with 1 sets of primers for methylated DNA and another set of primers for unmethylated DNA (Supplementary Table S4) designed by using MethPrimer website (https://www.urogene.org/methprimer/). The following PCR conditions were used for amplification: initial denaturation at 94°C for 5 min, followed by 40 cycles of 94°C for 30 s, annealing temperature of primers (53°C–58°C, varied among the primers) for 30 s and 72°C for 60 s; and final elongation for 7 min at 72°C (Ku et al., 2011). Amplified DNA products were electrophoresed for 25 min at 100V in a 2% Agarose gel. The gel was scanned, and the presence or absence of DMRs was determined based on the existence and strength of bands on the gel.

R package Methylkit, gplots, Dplyr, and Graph Pad Prism version 5.3.0 were used for PCA, Clustering, Heatmap generation and preparation of other figures. BAM files (*.bam) produced by Bismark were uploaded to SeqMonk software (http://www. bioinformatics. babraham.ac.uk/projects/seqmonk/) to visualize the location and methylation levels of genes. Data are reported as mean ± SEM. Two-tailed unpaired Student’s t-test was applied for the identification of significant difference (p < 0.05).

Semen straw samples produced during the late spring and winter from 5 different bulls were collected by ABS Global, of which sperm from three bulls in late spring and winter, respectively, were randomly selected and sequenced using WGBS (Figure 1A). The sequence statistics and alignment statistics are summarized in Table 1. WGBS generated an average of 112.3 ± 3.41 million reads (Table 1) with an average Phred quality score of 36% which indicates an overall high sequence quality. The bisulfite conversion rate averaged 99.7% ± 0.02%. The sequence alignment was conducted using bovine reference genome ARS-UCD 2.1 incorporated with Y chromosome Data (GenBank: CM011803.1), which exhibited 79.5%–82.7% of uniquely mapped reads across all the samples (Table 1). An average of 55.9 million reads were mapped uniquely with an average of 70% reads aligned on CpGs. Additionally, we detected on average 2.83% ± 0.01% of CHG and 3.24% ± 0.02% of CHH were methylated in sperm samples (Figure 1B). Sequence and alignment statistics listed in Table 1 did not show statistical differences when comparisons were made between winter and late spring groups. In brief, high quality sequencing data on DNA methylation were obtained from all sperm samples, which facilitated the subsequent downstream analysis.

For further analysis, we solely focused on methylation data on CpG context in order to obtain an overall scenario of genome-wide CpG methylation. In accordance with prior WGBS investigations of cow sperm samples, we found a genome-wide CpG methylation rate of 62.1% ± 2.3% and 60.7% ± 0.79% during the late spring and winter, respectively (Zhou et al., 2018), with no difference (p = 0.329) across two seasons (Figure 2A). To confirm, we performed an unsupervised hierarchical clustering analysis (Figure 2B) and principal component analysis (PCA) (Figure 2C). Clustering could not segregate samples according to the treatment groups and methylation difference among the individual samples was higher than the difference between groups. Additionally, PCA did not segregate the samples into two different groups, demonstrating that seasons had no impact on the global DNA methylation levels of sperm.

The CpGs covering at least ten uniquely mapped reads were filtered for further analysis by MethylKit and named CpG10 (Akalin et al., 2012). Around 30% of CpGs fell under the category of CpG10 and had an average methylation percentage of 26.3%. Among CpG10s, 33% had a methylation level <20% (hypomethylated), more than 60% had methylation level within 20%–80% (intermediate) and around 10% had methylation level more than 80% (hypermethylated) (Table 2). But the number of hypomethylated, intermediate and hypermethylated CpG10s did not differ significantly between late spring and winter groups.

Next, we ran unsupervised hierarchical clustering (Supplementary Figure S1A) and PCA (Supplementary Figure S1B) of filtered CpG10s. Consistent with previous findings, PCA and clustering could not distinguish samples into two different groups. Taken together, cold exposure of bulls during winter did not induce changes in overall methylation of CpG10s and the effect of inter-individual variability was higher than the effect of cold exposure.

To reveal the effect of seasons on DNA methylation, we analyzed 92,622 CpG10s obtained by filtering across 6 sperm samples. We applied a moderate threshold of DNA methylation difference ≥10% and p < 0.05 and R Bioconductor package Methylkit was able to detect 22 DMCs. Furthermore, we obtained 3,163 DMCs with a threshold of methylation difference ≥10% and q-value < 0.05 (Figure 3A; Supplementary Table S1) by R Bioconductor package DSS (Feng and Wu, 2019). Among the identified DMCs by DSS, methylation difference ranged from 12% to 82% across two groups and around 90% of the DMCs had methylation difference between 25% and 100% (Figure 3A). Therefore, we used a stricter threshold of q-value < 0.001 and methylation difference ≥25%, and only 8% of the previously identified DMCs fallen into this category. In the rest of the study, we used the threshold of ≥10% difference and q-value < 0.05 for identification of DMCs. The ratio of the hypermethylated and hypomethylated DMCs were 57.8% and 42.2% respectively (Figure 3B). Distribution of DMCs also varied between hypomethylation and hypermethylation groups, with ∼63% DMCs showed more than 50% difference in methylation percentage and ∼27% DMCs had an over 25% difference in the percentage of methylation (Figure 3B).

All the DMCs were then annotated with R Bioconductor package “Genomation” relative to the gene and gene features (Figure 3C). The majority of DMCs (93%) were spread over intergenic regions and introns while only 4% overlapped with the promoters. For CpG features, around 52% of DMCs were annotated outside of the CpG Islands (CGIs) or CpG shores (4 kb on both sides of CGI) while around 41% of the DMCs overlapped with CGIs. Furthermore, DMCs were spread over all the chromosomes (Supplementary Figure S2)

Differential methylation analysis was performed by using DSS. We used default DMR calling function “callDMR” of DSS to identify significant DMRs with a q-value < 0.05. DSS identifies several significant DMCs located close to each other and merge them to form a DMR. DSS was able to detect 438 DMRs (Supplementary Table S2), and the majority of the DMRs were hypomethylated (69.2%) (Figure 4A). The DNA methylation percentage of hypomethylated and hypermethylated DMRs was able to cluster samples according to two seasons on the heatmap (Figures 4D, E).

DMRs were then annotated to identify their overlapping with introns, exons, promoters, intergenic regions, CGIs, and shores. Among the hypermethylated DMRs, 60% overlapped with intergenic region and 13% overlapped with promoters. Regarding regions overlapped with CpG density, 34% overlapped with CGIs (Figure 4B). The annotation of hypomethylated regions exhibited similar trends (Figure 4C). Most hypomethylated DMRs overlapped with intergenic regions while only 7% overlapped with promoter regions. In terms of CpG density, most DMRs overlapped with CGIs (46%) (Figure 4B).

Hypo and hypermethylated DMRs were dispersed on almost every chromosome except chromosomes 26 and 28 which did not contain DMRs. Surprisingly, 50% of the hypomethylated DMRs were present on four chromosomes: 13, 16, 18 and Y, while Y chromosome contained 92 hypomethylated DMRs (Figure 5).

We performed the functional enrichment analysis in order to determine the biological processes and molecular functions of differentially methylated genes. We used Database for Annotation, Visualization and Integrated Discovery (DAVID) for functional enrichment analysis of 186 differentially methylated genes overlapped with 438 DMRs, among which 142 genes were distributed in 17 clusters among which top 3 clusters (Figures 6A–C) were associated with biological processes such as chemotaxis, transcription regulation from RNA polymerase II promoter, and histone modification.

We also did a Gene Ontology (GO) Enrichment Analysis on 186 genes under the category of biological processes, cellular components, and molecular functions. Several GO terms identified by GO enrichment analysis matched those identified by DAVID. Several DMGs were associated in crucial biological processes: RNA polymerase II transcription regulation, osteoblast differentiation, regulation of neuroblast proliferation, limbic system development, regulation of embryonic mitotic cell cycle and chemotaxis (Figure 7). Along with this, DMGs were enriched in top molecular functions: chemokine receptor binding, histone methyltransferase activity, histone reader activity and G-protein-coupled receptor binding (Figure 7); and cellular component terms: mitochondria, intercellular organelle lumen, cytoskeleton, and cytoplasm (Figure 7).

Reactome pathway analysis was performed to determine the pathways in which the DMGs were involved (Figure 7). “Post-translational protein modification”, “Histone Methylation”, “Signaling by Nuclear Receptors”, “ESR-mediated signaling”, “Transcriptional regulation by RUNX1” and “Wnt Signaling pathway” were among top 20 identified significant pathways.

DAVID and GO enrichment analysis both demonstrated that, a large portion of DMGs are involved in common GO terms: chemotaxis, transcription regulation, histone modification and embryonic osteogenesis, neurogenesis and mitotic cell cycle. To summarize, the Reactome pathway analysis and GO enrichment analysis revealed that methylated genes might play crucial roles in embryonic neuroblast and osteoblast differentiation through the Wnt Signaling pathway, transcriptional regulation, and histone modification.

We performed a literature review of 186 differentially methylated genes (DMGs) and identified eight unique DMGs: Pax6, Macf1, Mest, Ubqln1, Smg9, Ctnnb1, Lsm4, and Peg10 involved in embryonic development and nine unique DMGs: Prmt6, Nipal1, C21h15orf40, Slc37a3, Fam210a, Raly, Rgs3, Lmbr1, and Gan involved in osteogenesis in both animal and human studies (Table 3; Supplementary Figure S3).

Among the genes involved in embryonic development, Macf1, Ubqln1, Smg9, Ctnnb1 were hypomethylated and Pax6, Mest, Lsm4, and Peg10 were hypermethylated. The identified DMRs overlapped with promoters of Ubqln1, Smg9, and Peg10 and introns of rest genes (Table 3). Among the genes related to embryonic development, Mest and Peg10 are two paternally expressed imprinted genes that exhibited ∼50% increase in the methylation level.

Among the osteogenic genes, Prmt6, Nipal1, C21h15orf40, Slc37a3, and Lmbr1 were hypomethylated in winter whereas others were hypermethylated. The identified DMRs had 30%–55% difference in methylation profile across two seasons and overlapped with promoter regions of Prmt6, Nipal1, C21h15orf40, and Slc37a3; in other genes, DMRs overlapped with the CGI of exons in seven genes and introns of three genes (Table 3).

To verify the location of the DMRs identified by WGBS, we performed MS-PCR of three randomly selected genes known to regulate embryonic development: Smg9, Mest, and Peg10, and three genes associated with osteogenesis: Lmbr1, C21h15orf40, and Prmt6. The bisulfite converted DNA was amplified with methylated and unmethylated primers designed by targeting specific regions of DMRs and electrophoresed in Agarose gel. Differential methylation quantified in WGBS and bands found from MS-PCR of the above-mentioned genes are compared in Figures 8, 9. In WGBS, Peg10 and Mest had hypermethylated DMR with ∼30% higher methylation during winter which overlapped with their promoters’ region. In MS-PCR both genes showed stronger bands with winter samples amplified by methylated primers and strong bands in late spring samples with unmethylated primers (Figures 8A, C). On the other hand, Smg9, Lmbr1, C21h15orf40 and Prmt6 genes exhibited hypomethylated DMRs during winter. In MS-PCR all these genes had weak to no band during winter when amplified with methylated primers and stronger band during late spring amplified with unmethylated primers (Figure 8B; Figures 9A–C). The intensity of bands obtained from methylated and unmethylated primers corresponds with the results of WGBS and verifies the presence and absence of DMRS across two seasons. In conclusion, cold exposure induced differential methylation in genes associated with embryonic development, and osteogenesis and the identified DMRs overlapped with promoters and introns of genes including imprinted genes.