Differentially methylated CpG sites related to fertility in Japanese Black bull spermatozoa: epigenetic biomarker candidates to predict sire conception rate.

For semen suppliers, predicting the low fertility of service bull candidates before artificial insemination would help prevent economic loss; however, predicting bull fertility through in vitro assessment of semen is yet to be established. In the present study, we focused on the methylated CpG sites of sperm nuclear DNA and examined methylation levels to screen new biomarkers for predicting bull fertility. In frozen-thawed semen samples collected from Japanese Black bulls, for which the sire conception rate (SCR) was recorded, the methylation level of each CpG site was analyzed using human methylation microarray. According to regression analysis, 143 CpG sites related to SCR were significantly differentially methylated. Whole genome bisulfite sequence data were obtained from three semen samples and the differentially methylated regions (DMRs) that included the target CpG sites selected by human methylation microarray were confirmed. Using combined bisulfite restriction analysis, fertility-related methylation changes were detected in 10 DMRs. With the exception of one DMR, the methylation levels of these DMRs were significantly different between groups with high fertility (> 50%) and low fertility (< 40%). From multiple regression analysis of methylation levels and SCR, three DMRs were selected that could effectively predict bull fertility. We suggest that these fertility-related differences in spermatozoal methylation levels could be new epigenetic biomarkers for predicting bull fertility.

To protect the cattle breeding industry, the current problem of declining conception rates after artificial insemination (AI) requires urgent attention [1, 2]. For this purpose, i.e., to reduce risk, improve herd fertility, and minimize economic loss, cost-effective procedures such as the bull breeding soundness evaluation have been implemented. In Japan, beef cattle breeding focused on meat quality and using genetically biased bulls has seen an annual increase in the inbreeding coefficient of Japanese Black cattle, which has likely contributed to infertility issues [3]. In addition, summer heat stress due to the effects of global warming in Japan has led to summer sterility and reduced semen quality [4]. Indeed, many service bull candidates with excellent genetic traits are suffering from poor semen quality or even infertility. Therefore, predicting candidate bulls with low fertility before AI is desired by semen suppliers to prevent significant economic loss [1]. When attempting to predict bull fertility through in vitro semen assessment, sperm nucleus DNA damage and DNA methylation level abnormality, which are thought to be caused by environmental stress, should be considered as factors affecting low conception [5, 6]. DNA methyltransferase proteins are known to contribute to sperm DNA methylation and male fertility in general [7]. In humans, aberrant DNA methylation patterns have been correlated with abnormal semen parameters, idiopathic male infertility, and even pregnancy failure [6, 8,9,10]. In bulls, fertility status may also be associated with differences in sperm DNA methylation patterns [11]. In previous studies, bull sperm DNA methylation signatures were characterized by performing whole-genome DNA methylation binding domain sequencing [12], whole-genome bisulfite sequencing (WGBS) [13], and by detecting the variable methylated regions associated with the reproductive traits of bulls.

In our previous studies, we analyzed the genome-wide DNA methylation profiles of bull spermatozoa using a human DNA methylation microarray; accordingly, we identified the differential methylation regions (DMRs) with age-dependent different methylation (ADDM) levels [14,15,16]. Using combined bisulfite restriction analysis (COBRA), nine CpG sites were confirmed with nine ADDM-DMRs. These age-dependent methylation changes occurred rapidly at a young age and contributed to the integrity of spermatozoa and their fertility.

In the present study, we again analyzed the genome-wide DNA methylation profiles of bull spermatozoa using a human DNA methylation microarray and thereby determined the CpG sites with methylation levels in bulls classified into high or low fertility groups based on their sire conception rate (SCR). Methylation changes in some CpG sites can be visualized using COBRA [14, 15], and the methylation levels of these CpG sites can be analyzed from many samples. Using these techniques, we detected candidate fertility-associated DMRs (cFA-DMRs), for which the methylation levels were related to SCR, to screen for new biomarkers that predict bull fertility.

Materials and Methods

All chemicals used in the study were purchased from Sigma-Aldrich (St. Louis, MO, USA) unless otherwise indicated.

Semen sampling

Cryopreserved semen samples (n = 50) were collected from 28 Japanese Black bulls kept at four facilities (Hiroshima Prefectural Livestock Technology Research Center, Gifu Prefectural Livestock Research Institute, Tottori Prefectural Livestock Research Institute, and a private ranch) located in different areas of Japan. These samples were used to investigate conception rates after AI (i.e., SCR, %). The ages of the bulls from which the samples were taken and used for DNA methylation analysis are shown in Fig. 1. To analyze the age effect, cryopreserved semen samples (n = 73) in the age range 10–162 months were collected from 38 bulls kept at five facilities (adding Ibaraki Prefectural Livestock Research Center) and then used for subsequent COBRA.

Fig. 1.

Sire conception rate (SCR, %) and age (in months) of each semen sample used for EPIC microarray analysis (A: n = 17 from 14 bulls) and combined bisulfite restriction analysis (COBRA; B: n = 50 from 28 bulls). Samples were divided into high fertility (HF: > 50% fertility) and low fertility (LF: < 40% fertility) groups. The LF samples categorized in the LF1 cluster by EPIC (see Fig. 2C) are shown as white squares; the LF sample that was categorized in the HF cluster by EPIC (see Fig. 2C) is shown as an asterisk.

Fig. 2.

Heatmaps (A and B) showing the different methylation levels analyzed by EPIC microarray at (A) 332 selected CpG sites that differed between high fertility (HF, n = 6) and low fertility (LF, n = 7) sample groups and (B) 143 differential CpG sites selected according to regression analysis using a linear model of each CpG site in 17 samples with sire conception rate data. The dendrograms above and to the left of the heatmaps show hierarchical clustering results among individual animals and among the CpG sites based on the methylation data, respectively. The gray scale at the upper left of each heatmap indicates methylation status, with white (0.0) being the lowest and black (1.0) being the highest methylation level. A cluster analysis (C) was also conducted using the methylation levels of the 143 CpGs shown in B. The LF sample categorized into the HF cluster is shown by an asterisk.

DNA methylation analysis using the Infinium EPIC BeadChip array

Thawed semen samples (n = 17 from 14 bulls; Fig. 1A) were washed twice in Dulbecco’s phosphate-buffered saline without calcium chloride or magnesium chloride. Genomic DNA was extracted from the semen samples using a DNeasy Blood and Tissue Kit (Qiagen, Valencia, CA, USA) according to the manufacturer’s instructions but with some modifications as previously described [14].

The Infinium MethylationEPIC BeadChip (Illumina, Inc., San Diego, CA, USA) array was used for genome-wide bovine DNA methylation analysis as previously described [15]. Briefly, the processing of DNA for the methylation arrays was performed according to the Illumina protocols at Takara Bio Inc. (Shiga, Japan). After methylation array analysis, informative CpG sites were selected, and the fluorescence signal intensities of the methylated and unmethylated alleles were obtained [17]. Three analyses were performed as follows with the R software (https://www.R-project.org/). In Experiment 1, the different methylation levels of the CpG sites were selected by comparing results between samples with high fertility (HF; > 50% fertility; n = 3, H1–3) and low fertility (LF; < 40% fertility; n = 3, L1–3) using a Welch’s two sample t-test for each CpG site. In Experiment 2, the CpG sites were similarly selected by comparing HF (n = 6; H1–4, H5-1, and H5-2) and LF (n = 7; L1–4, L5-1, L5-2, and L6) samples using a Welch’s two sample t-test and regression analysis of the methylation rate with respect to the conception rate for each CpG site. In Experiment 3, the CpG sites were selected by regression analysis using a linear model of each CpG site in 17 samples (H1–4, H5-1, H5-2, H6, H7-1, H7-2, L1–4, L5-1, L5-2, L6, and L7) with SCR data. The most reliable CpG sites were quality checked and those with detection P values = 0 (defined by Illumina) and total signal intensities > 1,000 were selected for each experiment. Ultimately, candidate CpG sites with different methylation levels associated with fertility were identified.

To select suitable CpG sites in the DMRs, WGBS analysis was performed for three samples (L5-1, L8, and H7-1; the conception rate after AI was 35.7%, 37.8%, and 62.2%, while age was 13, 14, and 14 months, respectively), and neighboring CpG sites in DMRs were defined using the results of WGBS. Each methylation level of the CpG site located in the range 2,000 base pairs (bp) upstream and downstream of the target CpG site selected by microarray analysis was then aligned. Subsequently, the candidate CpG sites with different methylation levels associated with bull fertility were listed and compared among individuals. From these candidate sites, suitable CpG site restriction enzymes for methylation level determination were selected for use in subsequent COBRA.

Detection of differentially methylated regions by COBRA

COBRA was performed according to a previously published method [15]. Briefly, bovine genome sequences (ARS-UCD1.2) were identified using BLAST (http://blast.ncbi.nlm.nih.gov/Blast.cgi) from human probe sequences and then bisulfite primers were designed using MethPrimer (http://www.urogene.org/methprimer/index.html). Genomic DNA extracted from frozen-thawed semen samples was bisulfite-converted using a MethylEasy Xceed Rapid DNA Bisulphite Modification Kit (Human Genetic Signatures Pty. Ltd., New South Wales, Australia) and then amplified by polymerase chain reaction (PCR) using a TaKaRa EpiTaq HS kit (Takara Bio Inc.) with the following conditions: 94°C for 2 min followed by 35 cycles of 94°C for 30 sec, 50–58°C (depending on the primer sets; Table 1) for 30 sec, and 72°C for 30 sec. PCR products were digested with the restriction enzymes shown in Table 1. The digested fragments were analyzed using 3% agarose gel electrophoresis. The intensities of the digested and undigested fractions were measured by densitometry using ImageJ 1.53a (http://imagej.nih.gov/ij; National Institutes of Health, Bethesda, MD, USA), and methylation levels were calculated based on the ratio between the intensities of the digested and undigested DNA band. The DMRs associated with fertility were selected for subsequent analysis.

Table 1.

Primers and locations of target CpGs revealed by combined bisulfite restriction analysis (COBRA)

ID	Forward and reverse primer sequences	Product length (bp)	Annealing (°C)	Restriction enzyme	No. of restriction sites	Restriction fragment length a (bp)	CpG location (chromosome no.)	Features	Target ID of EPIC
CpG-F1	5'TAAATGGTTTTAGTAAGAAATTAAATATAA3'	196	50	HpyCH4IV	2	160, 28, 8	NC_037353.1, 38105097 (26)	5': empty spiracles homeobox 2	cg13647079
	5'CCAACTAAACAAATCATTATAAAACTA3'							3': RAB11 family interacting protein 2
CpG-F2	5'GGTTTTGTGTGGTTGTATAGTGAAT3'	188	50	BstUI	1	141, 47	NC_037334.1, 11271815 (7)	DnaJ heat shock protein family (Hsp40) member B1	cg07483523
	5'AAATCCTTTCAAAAAAAATCACTTC3'
CpG-F3	5'GGGTATTGGGGAGATATTTTTGTAT3'	180	54	TaqI	2	110, 35, 35	NC_037350.1, 27193986 (23)	Neurogenic locus notch homolog protein 4 isoform X1neurogenic locus notch homolog protein 4 precursor	cg08801479
	5'CAACCCCAATATACACTAACCTAACA3'
CpG-F4	5'GAAAATGTGGGAAAATTTATATTTTTG3'	283	56	BstUI	1	165, 118	NC_037350.1, 22717962 (23)	-	cg23797553
	5'ATCAACCATCCATCCATCTAATTAA3'
CpG-F5	5'TAGAGAGGTTATTTGGGAGGTATTT3'	215	52	TaqI	1	139,76	NC_037348.1, 65954056 (21)	5' side: retrotransposon-like protein 1	cg05537796
	5'CCCTTCAAATACCAAAAAAAATACTAA3'							3' side: thyroxine 5-deiodinase
CpG-F6	5'TGTTTTTTTGGTATGGTTTTTTGTT3'	228	50	AciI	1	151, 77	NC_037343.1, 4869926 (16)	Ubiquitin thioesterase OTU1	cg07483523
	5'AACATCCACTATAATCCACTTCATCTTAT3'
CpG-F7	5'GGGGGTTTAGTTTTTTAGTTTTTTAAAT3'	136	56	TaqI	1	100, 36	NC_037353.1, 24169514 (26)	SH3 and PX domain-containing protein 2A isoform X2SH3 and PX domain-containing protein 2A isoform X3	cg04762698
	5'CACCACATACAATACCTACCAAAAA3'
CpG-F8	5'GTGGAGTTTGGGTTATTTATTTTTG3'	213	57	TaqI	1	162, 51	NC_037350.1, 7490261 (23)	Death domain associated protein	cg27584448
	5'CACCACCTCTAATAAACCCTCTAAA3'
CpG-F9	5'TTTTGTAGGTAGGAAGTTGGATTGT3'	290	52	HpyCH4IV	1	227, 63	NC_037338.1, 10665097 (11)	Methylcytosine dioxygenase TET3 isoform X1	cg17184593
	5'CAACTAATACAAACCCCACAATTTT3'
CpG-F10	5'TGAATAGTAGTTGAATATGGGTTAGT3'	236	50	AciI	3	134, 102	NC_037333.1, 109018628 (6)	5' side: biorientation of chromosomes in cell division 1 3' side: cytoplasmic polyadenylation element binding protein 2	cg05696584
	5'CCCTACCCTTCACAAAAAAAA3'

a) Restriction fragment length when all CpGs located in restriction endonuclease cleavage sites are methylated.

Analysis of methylated levels in cryopreserved semen by COBRA

DNA extracted from cryopreserved semen samples (n = 50 from 28 bulls; SCR 0%–70%; age: 12–144 months; Fig. 1B) was used for COBRA. The methylation levels of the selected DMRs were estimated by COBRA and their correlation with fertility was analyzed for each CpG by Spearman’s rank correlation using GraphPad Prism (for Windows, version 7). The differences in methylation levels between the HF and LF groups (defined by SCR: HF > 50%, n = 21; LF < 40%, n = 20) were analyzed by Mann–Whitney U test using GraphPad Prism.

Prediction of bull fertility using the methylation levels analyzed by COBRA

Bull fertility (i.e., estimated SCR, %) was evaluated using the methylation levels of the 10 CpG sites in the sperm DNA analyzed by COBRA. First, multiple regression analysis for each methylation level of all 10 CpG sites was performed using the SCR as an objective variable and methylation levels at the 10 CpG sites as explanatory variables. Subsequently, a stepwise forward and/or backward method based on Akaike’s Information Criterion (AIC) was performed to select effective CpG sites to obtain a simple and easily interpretable model. A predictive formula for estimating bull fertility was then created using selected CpG sites. Finally, the estimated SCR was calculated from the methylation levels of all or selected CpG sites, and this value was compared with the actual SCR. In addition, as age was related to bull fertility, variable selection was performed that included age in months as a variable.

Results

DNA methylation status revealed by the Infinium MethylationEPIC BeadChip array

In Experiment 1, 84,954 CpG sites passed the data quality control process in six examined samples (9.81%). Of these sites, 349 showed significant differences in beta values between the HF and LF groups (beta-value difference cutoff = 0.1). In Experiment 2, 79,003 CpG sites passed the data quality control process in all examined samples (9.12%). Of these, 332 CpG sites showed significant differences in beta values between the analyzed groups (beta-value difference cutoff = 0.1; Fig. 2A). In Experiment 3, 73,442 CpG sites passed the data quality control process in all examined samples (8.48%). Using a cutoff value of 0.3 for differences in beta values between the “MAX” and “MIN” values analyzed (n = 17), 143 CpG sites showed significantly different beta values (Fig. 2B). The correlation coefficient between the average estimated SCR, according to the regression equation from the methylation rate of these CpG sites, and the actual SCR was r = 0.87. In a cluster analysis of the methylation levels from the 143 CpGs, part of the LF group branched (LF1) before the LF and HF groups branched (LF2 and HF), i.e., three clusters were produced (Fig. 2C). However, one of the LF samples was divided into a HF cluster (Figs. 1 and 2). The LF1 cluster included samples with a semen collection age ≤ 20 months, as shown in Fig. 1.

Candidate DMRs associated with bull fertility confirmed by COBRA

In Experiment 1, 57 matching bovine DNA sequences from 349 differentially methylated CpG sites (16.3%) were selected after BLAST analysis. Thirteen CpG sites, including restriction enzyme sites, were selected and primer sets were designed. Three of these 13 sites were confirmed as DMRs by COBRA; one site (CpG-F1, Table 1) was identified as a cFA-DMR from further screening. In Experiment 2, 173 matching bovine DNA sequences were obtained from 332 differentially methylated CpG sites (52.1%) following BLAST analysis. COBRA primer sets were designed for each of the 14 selected CpG sites. Four of these were confirmed as DMRs by COBRA and three (CpG-F2–4, Table 1) were identified as cFA-DMRs. In Experiment 3, 48 of 143 probe sequences including CpG sites (33.6%) were matched with bovine DNA sequences after BLAST analysis. Twelve CpG sites, including restriction enzyme sites, were selected and primer sets were designed. Two of these 12 sites were also detected in Experiment 2 (CpG-F2 and CpG-F3). From the remaining 10 regions, four were confirmed as DMRs by COBRA and two were identified as cFA-DMRs (CpG-F5 and CpG-F6, Table 1).

Many CpG sites associated with fertility according to EPIC analysis were inapplicable to the COBRA method. DMRs including the target CpGs were defined by comparison with peripheral sequences (4 kbps) analyzed by WGBS (Fig. 3); thus, CpGs that showed single nucleotide polymorphisms (SNPs) were excluded. Therefore, four additional CpGs located near the CpGs detected by EPIC were used for COBRA analysis; the regions including these CpGs were identified as cFA-DMRs (CpG-F7–F10, Table 1).

Fig. 3.

The methylation differences of peripheral regions (4 kbps) containing the target CpG sites revealed by whole genome bisulfite sequencing results (n = 3; L5-1, L7, and H7-1). The peripheral areas of CpG-F1, -F3, -F4, -F5, -F8, and -F9 are shown. Arrows indicate target CpG sites identified by EPIC microarray analysis.

Fertility-related methylation changes in the 10 identified cFA-DMRs were confirmed by COBRA using independent samples (n = 50, Table 2 and Fig. 4). As shown in Fig. 4, a regression line generated for each cFA-DMR indicated that methylation levels at the cFA-DMRs increased as SCR increased (Fig. 4). Furthermore, there was a significant difference in methylation levels between the HF and LF groups at each cFA-DMR (Table 2, see “HF vs. LF”), with the exception of CpG-F6 (P = 0.051). In addition, there was a significant correlation between the methylation level at each cFA-DMR and the age of the sample (P < 0.01). In linear regression analysis incorporating the age effect, the methylation levels of eight cFA-DMRs (namely CpG-F1–F5 and CpG-F8–F10) remained significantly associated with SCR (P < 0.05; Table 2, see “SCR-Age”).

Table 2.

Effectiveness of candidate fertility-associated–differential methylation regions (cFA-DMRs) revealed by significant differences (P < 0.05) between methylation levels and sire conception rate (SCR, %) or between fertility groups [high fertility (HF) vs. low fertility (LF)]

	EPIC		COBRA
ID	SCR a	HF vs. LF c	SCR a	HF vs. LF c	Age a	SCR-Age a,b
CpG-F1	0.578	0.447	< 0.001	< 0.001	< 0.001	< 0.001
CpG-F2	0.007	0.013	< 0.001	< 0.001	< 0.001	0.001
CpG-F3	0.018	0.010	< 0.001	< 0.001	0.036	< 0.001
CpG-F4	0.054	0.103	0.004	0.031	< 0.001	0.040
CpG-F5	0.002	0.009	< 0.001	0.008	< 0.001	0.002
CpG-F6	0.005	0.043	0.017	0.051	< 0.001	0.195
CpG-F7	0.080	0.118	0.015	0.028	0.012	0.102
CpG-F8	0.013	0.016	< 0.001	0.002	< 0.001	0.002
CpG-F9	0.008	0.004	0.002	< 0.001	< 0.001	0.036
CpG-F10	0.005	0.015	0.002	< 0.001	0.004	0.023

Results of EPIC microarray analysis and combined bisulfite restriction analysis (COBRA) are shown. a) For SCR, Age, and SCR-Age, differences were analyzed by linear regression. b) For SCR-Age, linear regression analysis was performed while considering the effect of age. c) For HF vs. LF, differences between the HF and LF groups were analyzed using a Mann-Whitney U test.

Fig. 4.

Fertility-related methylation changes identified as candidate fertility-associated–differential methylation regions (cFA-DMRs; i.e., CpG-F3–F5) according to combined bisulfite restriction analysis (COBRA). Upper panels: according to linear regression analysis, the correlations between the methylation level of each DMR and the sire conception rate (SCR, %) of the sample were statistically significant (P < 0.01, n = 50). The equation and R2 value from the analyses have been reported with each graph. Lower panels: significant differences were observed between the methylation levels of the high fertility (HF, n = 21) and low fertility (LF, n = 20) groups (Mann-Whitney U test, * P < 0.05, ** P < 0.01).

From stepwise selection based on AIC, three of the cFA-DMRs, namely CpG-F3 (P < 0.01), CpG-F4 (P < 0.05), and CpG-F5 (P = 0.05), were selected as the main effective contributors (r = 0.67) to the prediction of SCR (in comparison to all 10 cFA-DMRs for which r = 0.71). SCR predictions were performed using the two multiple regression equation models and these revealed by the methylation levels of all 10 cFA-DMRs or these 3 selected cFA-DMRs (Fig. 5). Correlations between the actual conception rate and the predicted SCRs obtained using the methylation levels of 10 or 3 CpG sites yielded coefficients of 0.71 or 0.67, respectively (Fig. 5).

Fig. 5.

The predicted sire conception rate (SCR, %) for 50 samples calculated using multiple regression equations with the methylation levels of 10 candidate fertility-associated–differential methylation regions (cFA-DMRs; left panel) and three selected cFA-DMRs (CpG-F3–F5; right panel). Correlations between the actual and predicted SCRs produced coefficients of 0.71 and 0.67 for the 10 and 3 cFA-DMR analyses, respectively. The equation and R2 value from the analyses are reported with each graph.

Discussion

Although bull breeding soundness examinations and semen analyses can be used to predict bull fertility, an effective and low-cost screening method is also required to exclude LF bulls from breeding programs [1]. Indeed, for this purpose, it is not sufficient to simply evaluate physical soundness and semen quality, i.e., via qualities such as progressive motility, morphological normality, and acrosome integrity. Bull fertility can usually be confirmed by analysis of SCR in cows after AI, and SCR is also an important factor in the selection of service bulls; however, it is difficult to incorporate SCR into fertility prediction when selecting young bulls as service bull candidates. Although semen property and in vitro fertilization (IVF) tests are effective, they are not sufficient to identify sires with low SCRs [12]. In recent years, protein biomarkers used to detect subfertile bulls have been reported, e.g., sperm acrosome-associated 1 proteins [18, 19]. The relationship between epigenetic data and bull fertility has also been investigated to identify new epigenetic biomarkers [11, 12, 20]. These studies, however, have been difficult to apply to Japanese Black bulls, perhaps because of the different methylation statuses between bovine breeds with different genomic backgrounds that were raised in different countries [21]. Indeed, significant differences in the DNA methylation levels of sperm in bovine breeds have been observed among breeds on a global scale [21]. It is possible that the presence of DNA polymorphism affects global CCGG content and methylation.

In the present study, we screened biomarkers for SCR prediction using the methylation levels at CpGs detected by the human platform EPIC and we applied our findings to the bovine genome. Specifically, we identified 10 cFA-DMRs in Japanese Black bulls using the COBRA method. LF bulls were identified using 332 CpGs with methylation levels that differed between the HF and LF groups (according to EPIC detection) and via cluster analysis of these groups. Thus, microarray analysis of semen samples can apparently reflect bull fertility (i.e., HF and LF groups) in bulls for which fertility is unknown.

Using the 143 CpG sites shown to have a high correlation between SCR and methylation level, along with cluster analysis, the HF and LF groups were almost completely separated. The LF group was divided into two clusters, one close to the HF group (LF2) and more distant from this group (LF1, which contained semen of ≤ 20 months of age). One individual was divided into LF1 (L5-1, ≤ 20 months of age) and LF2 (L5-2, > 20 months of age), which may reflect differences in age. Thus, our results reveal that the methylation levels of CpGs change according to the bull’s age. Therefore, when using ADDM-DMRs as an evaluation index for SCR, it is necessary to consider age-related changes in sperm DNA methylation status.

We found nine ADDM-DMRs that had previously been reported [15]; four of these were associated with SCR (data not shown) and were not detected as cFA-DMRs in this study. We observed rapid changes in methylation levels until the bulls reached four years of age; therefore, the methylation level of bull sperm DNA is apparently related to sexual maturation. A previous study showed that spermatozoa from bulls aged 10 months had a unique epigenetic profile that could compromise their value [22]. In another study, varying nutrition in prepubertal calves resulted in some changes to their postpubertal sperm DNA methylation profile [23]. The age at which semen collection is possible and at which semen properties are sufficiently stable for use in breeding differs depending on the bulls. Therefore, when predicting fertility in bulls, it is important to consider not only the methylation level and age but also semen characteristics (which can be evaluated using a standard method).

Having used the methylation levels of 10 cFA-DMRs measured by COBRA in a multiple regression equation for SCR prediction, the correlation coefficient between the actual and predicted SCRs was 0.71. Although the target CpGs could be increased by further EPIC analysis, only three of the 10 CpGs were effective and necessary for SCR prediction. Stepwise selection/regression is a way of selecting important variables for obtaining a simple and easily interpretable model. The rest seven CpGs did not significantly affect the prediction results (P > 0.05). This is because similar patterns of methylation changes were observed in 10 CpGs, e.g., the higher the conception rate, the higher the methylation level. Moreover, with age in months also being considered in the analysis, the methylation levels of the three effective cFA-DMRs (i.e., CpG-F3–F5) were still significantly associated with SCR. Five of the 10 cFA-DMRs (F1, F2, F5, F8, and F9) showed significant differences among facilities (data not shown). There was no significant difference of SCR among the facilities; however, there was a significant difference of age, indicating that the location of facilities may be reflective of the age effect. Therefore, these fertility-related differences in spermatozoa methylation levels could be new epigenetic biomarkers for bull fertility prediction. It should be noted that the multiple regression equation used for predictions must be constantly updated since the DMRs and coefficients selected by multiple regression can change depending on the specificity and degree of contribution according to the samples analyzed. This implies that it is important to continue updating valid epigenetic information, including information on whether these DMRs can be used with other Japanese Black bulls. Further study will be needed to analyze the effect of animal origin and the housing environments.

WGBS data confirmed that each detected CpG was located in a DMR, respectively. Some were shown as SNPs, so these CpGs were removed from analysis. Unfortunately, from the WGBS data, it was not possible to detect differences between HF and LF samples because our sample numbers were limited. However, WGBS is considered a more effective method for finding unknown DMRs than EPIC. To detect cFA-DMRs by WGBS, many more samples should be analyzed, which would be expensive. In contrast, evaluating samples using COBRA with the identified cFA-DMRs is a simple and inexpensive method that provides rapid results.

Capra et al. [24] reported that methylation varies between high and low motile bovine sperm populations, and that maintenance of chromosome structure through epigenetic regulation is probably crucial for ensuring sperm functionality. The findings suggest that the DMRs may reflect the proportion of several types of sperm present in a single ejaculation. Therefore, the effects of different methylation statuses in sperm DNA on embryonic gene expression should be studied in future research to explore the underlying mechanism of bull infertility.

In conclusion, our study reveals that the methylation levels at some DMRs observed in sperm DNA are related to SCR. We suggest that these cFA-DMRs could serve as new epigenetic biomarkers for predicting SCR in bulls.