Genetic and Genomic Analyses of Service Sire Effect on Female Reproductive Traits in Holstein Cattle

Fertility and reproductive performance are key drivers of dairy farm profitability. Hence, reproduction traits have been included in a large majority of worldwide dairy cattle selection indexes. The reproductive traits are lowly heritable but can be improved through direct genetic selection. However, most scientific studies and dairy cattle breeding programs have focused solely on the genetic effects of the dam (GED) on reproductive performance and, therefore, ignored the contribution of the service sire in the phenotypic outcomes. This study aimed to investigate the service sire effects on female reproductive traits in Holstein cattle from a genomic perspective. Genetic parameter estimation and genome-wide association studies (GWAS) were performed for the genetic effect of service sire (GESS) on conception rate (CR), 56-day non-return rate (NRR56), calving ease (CE), stillbirth (SB), and gestation length (GL). Our findings indicate that the additive genetic effects of both sire and dam contribute to the phenotypic variance of reproductive traits measured in females (0.0196 vs. 0.0109, 0.0237 vs. 0.0133, 0.0040 vs. 0.0289, 0.0782 vs. 0.0083, and 0.1024 vs. 0.1020 for GESS and GED heritability estimates for CR, NRR56, CE, SB, and GL, respectively), and these two genetic effects are positively correlated for SB (0.1394) and GL (0.7871). Interestingly, the breeding values for GESS on insemination success traits (CR and NRR56) are unfavorably and significantly correlated with some production, health, and type breeding values (ranging from −0.449 to 0.274), while the GESS values on calving traits (CE, SB, and GL) are usually favorably associated with those traits (ranging from −0.493 to 0.313). One hundred sixty-two significant single-nucleotide polymorphisms (SNPs) and their surrounding protein-coding genes were identified as significantly associated with GESS and GED, respectively. Six genes overlapped between GESS and GED for calving traits and 10 genes overlapped between GESS for success traits and calving traits. Our findings indicate the importance of considering the GESS when genetically evaluating the female reproductive traits in Holstein cattle.

Fertility and reproductive performance are key drivers of dairy farm profitability. Hence, reproduction traits have been included in a large majority of worldwide dairy cattle selection indexes. The reproductive traits are lowly heritable but can be improved through direct genetic selection. However, most scientific studies and dairy cattle breeding programs have focused solely on the genetic effects of the dam (GED) on reproductive performance and, therefore, ignored the contribution of the service sire in the phenotypic outcomes. This study aimed to investigate the service sire effects on female reproductive traits in Holstein cattle from a genomic perspective. Genetic parameter estimation and genome-wide association studies (GWAS) were performed for the genetic effect of service sire (GESS) on conception rate (CR), 56-day non-return rate (NRR56), calving ease (CE), stillbirth (SB), and gestation length (GL). Our findings indicate that the additive genetic effects of both sire and dam contribute to the phenotypic variance of reproductive traits measured in females (0.0196 vs. 0.0109, 0.0237 vs. 0.0133, 0.0040 vs. 0.0289, 0.0782 vs. 0.0083, and 0.1024 vs. 0.1020 for GESS and GED heritability estimates for CR, NRR56, CE, SB, and GL, respectively), and these two genetic effects are positively correlated for SB (0.1394) and GL (0.7871). Interestingly, the breeding values for GESS on insemination success traits (CR and NRR56) are unfavorably and significantly correlated with some production, health, and type breeding values (ranging from −0.449 to 0.274), while the GESS values on calving traits (CE, SB, and GL) are usually favorably associated with those traits (ranging from −0.493 to 0.313). One hundred sixty-two significant single-nucleotide polymorphisms (SNPs) and their surrounding protein-coding genes were identified as significantly associated with GESS and GED, respectively. Six genes overlapped between GESS and GED for calving traits and 10 genes overlapped between GESS for success traits and calving traits. Our findings indicate the importance of considering the GESS when genetically evaluating the female reproductive traits in Holstein cattle.

INTRODUCTION
In recent decades, the genetic selection for functional traits, including reproductive performance, has received great emphasis in dairy cattle selection indexes, aiming to achieve more balanced and sustainable breeding goals (Egger-Danner et al., 2015). Female fertility has only been broadly included in national selection indexes of dairy cattle breeding programs over the past five decades . Although most dairy cattle breeding programs focus only on the genetic effects of the dam (GED) when genetically evaluating fertility and calving traits, there are evidence that the service sire may also have a genetic influence on female reproductive performance (Barton et al., 1984;Van Tassell et al., 2003;Averill et al., 2004)-for instance, the direct effects of the service sire (e.g., semen quality and viability) on female reproductive performance have been considered as indirect effects (Jansen, 1985). Jaton et al. (2017) reported that the heritability estimates of service sire on embryo quality were lower than the donor (0.02 versus 0.04) but still statistically significant. In 2008, a national evaluation model of sire conception rate (SCR) was established in the United Sates by the Animal Improvement Program Laboratory of the United States Department of Agriculture (Norman et al., 2008). SCR is measured as confirmed pregnancy ratio (in percentage) of each service sire. They also implemented a sire-maternal grandsire (S-MGS) model to estimate the genetic component of service sire in calving performance Jiang et al., 2018). Even though low heritability estimates have been reported for indirect indicators of male fertility (Berry et al., 2011;Tiezzi et al., 2013), it can still be improved through genomic selection (Lillehammer et al., 2011). Hence, understanding the genetic mechanisms underlying male fertility and developing accurate genomic prediction models are of great importance but still underexplored.
The determination of the genetic effect of service sire (GESS) on reproductive traits relies on genetic and genomic analyses, including the estimation of genetic parameters and genome-wide association studies (GWAS). Previous studies have identified interesting genomic regions associated with SCR, as reviewed by Taylor et al. (2018). However, there are few reports of the GESS on other reproductive traits at the genomic level (Fang et al., 2019;Jiang et al., 2018). In this context, the main objectives of this study were as follows: (Egger-Danner et al., 2015) to investigate the genetic background of service sire on female reproductive performance, including conception rate (CR), 56-day non-return rate (NRR56), calving ease (CE), stillbirth (SB), and gestation length (GL) in Holstein cattle and  to identify genomic regions and candidate genes associated with GESS on female reproduction performance.

Phenotypic and Genomic Data
Phenotypes Records of birth dates, insemination events, pregnancy diagnoses, and calving information collected in 39 farms (Sunlon Livestock Development Co. Ltd., Beijing, China) from 1987 to 2020 were extracted from the AfiFarm software (AfiFarm 1 ) and used in this study. The reproductive traits include CR (1 = pregnant, 0 = non-pregnant), NRR56 (1 = non-return, 0 = return), CE (scores from 1 to 3, in which 1 = unassisted, 2 = easy pull, 3 = hard pull and surgery needed), SB (1 = calf was alive 24 h after birth and 2 = calf was dead), and GL (in days). The last insemination before a positive pregnancy diagnosis in each parity was considered as pregnancy. Furthermore, the cows with calving records but without a positive pregnancy diagnosis were considered pregnant when last inseminated before calving. Insemination records that had the next insemination within 1-17 days or had neither insemination nor calving records after that time period were excluded from the NRR56 calculation. CE and SB were directly coded from raw data after excluding ambiguous records (3%) caused by mis-recording. GL records lower than 260 and greater than 302 days were also deleted. A descriptive summary for each trait after data editing is shown in Table 1. In total, 163,818 Holstein cows had phenotypes and were serviced by 1,952 bulls. The average (±SD) number of services per bull was 489 ± 1,947. A pedigree file spanning over 13 generations was provided by the Beijing Dairy Cattle Center (BDCC, Beijing, China) and consisted of 503,118 cows and 151,273 bulls born between 1957 and 2020. The estimated breeding values (EBVs) for six production traits (milk yield, milk protein yield, milk protein ratio, milk fat yield, milk fat ratio, and somatic cell score), two health traits (reproductive diseases and udder diseases), and five type traits (dairy character, milking system, capacity, rump, and overall conformation) traits were also provided by BDCC.

Genotypes
A total of 3,477 Holstein bulls were genotyped with the Illumina BovineSNP50 BeadChip (Illumina, Inc., San Diego, CA, United States) containing 54,609 single-nucleotide polymorphism (SNP) markers. These genotypes were imputed to the Illumina 150K Bovine Beadchip (containing 123,268 SNPs) using the BEAGLE v5.1 software (Browning et al., 2018) and a reference population consisting of 3,119 cows and 81 bulls. The SNP information was updated from an older version of the cattle reference genome (UMD 3.1) assembly to the current one (ARS-UCD 1.2) using the UCSC LiftOver tool 2 . Eight thousand five SNPs with missing position in the latest reference genome were removed from further analyses. The reference population was divided into a subreference population (2,725 individuals genotyped with the CR, conception rate (0 or 1); NRR56, 56-day non-return rate (0 or 1); CE, calving ease (1, 2, or 3); SB, stillbirth (0 or 1), GL, gestation length (day); SD, standard deviation; N, number of records.
150K SNP panel) and a sub-validation population (475 individuals genotyped with the 150K SNP panel but masked to only the 50K panel SNPs) for assessing the accuracy of genotype imputation as the concordance rate of imputed SNPs. Only SNPs with imputation accuracy greater than 90% were kept for further analyses. Furthermore, SNPs with minor allele frequency lower than 0.05, unknown chromosome or genome position, and extreme deviation from the Hardy-Weinberg equilibrium (p-value lower than 10 −6 ) were removed. After data editing, 109,274 SNPs located in the autosomes and pseudo-autosomal regions of the X chromosome were retained in the dataset.

Estimation of Genetic Parameters
The variance and covariance components for each reproductive trait were estimated using the AI-REML algorithm implemented in the DMUAI module of the DMU v6 software (Madsen et al., 2006). A previous study evaluating similar traits indicated that the heritability estimates of linear and threshold models tend to be similar (Meijering, 1985;Boichard and Manfredi, 1994). Therefore, the following linear mixed model was fitted: where y represents the vector of phenotypic observations (i.e., CR, NRR56, CE, SB, or GL), β is the vector of fixed effects included in the model, in which different systematic effects were fitted for each trait [i.e., AI technician, parity, semen type, and number of inseminations for CR and NRR56 and calf sex, parity, and group of calf size (divided according to their birth weights: 30-40, 40-50, and 50-60 kg) for CE, SB, and GL]. All of the fixed effects significantly (P < 0.05) influenced the dependent variables (CR, NRR56, CE, SB, and SB) and were identified based on mixed model analysis using the PROC MIXED function implemented in the SAS software (version 9.1.3; SAS Institute Inc., Cary, NC, United States); u m and u f are the vectors of the random animal effects accounted by GESS and GED, respectively; pe m and pe f are the vectors of the random permanent environmental effects of the service sire and dam, respectively; hym is the vector of the random herd-yearmonth effects; e is the vector of the random residual effects; and, X, Z 1 , W 1 , W 2 and Z 2 are the corresponding incidence matrices. We assumed that: where σ 2 m and σ 2 f are the additive genetic variances of service sire and dam, respectively; σ m,f is the genetic covariance of service sire and dam; σ 2 pem and σ 2 pef are the permanent environmental variances of service sire and dam, respectively; σ 2 hym is the herd-year-month variance; σ 2 e is the residual variance; ⊗ is the Kronecker product function; A is the additive genetic relationship matrix among the animals; and I is an identity matrix. For calving traits, the model could be considered as an improved siredam model that assumes that the GESS has a genetic covariance with the GED. Compared to the traditional evaluation models of calving traits that only consider the animal effect, the current model includes both service sire and dam effects in the female reproductive performance phenotypes through GESS and GED, respectively. Differently from direct and maternal (paternal) effects usually assumed to be correlated, this was not the case for the genetic components of service sire and dam in the current study. Therefore, the GESS heritability estimates were calculated as follows: and the repeatability estimates were calculated as follows: The GED repeatability was estimated in the same way but replacing σ 2 m and σ 2 pem in the numerator by σ 2 f + σ 2 pef . The standard error of the heritability and repeatability estimates, respectively, were calculated using the Delta method (Su et al., 2007). A Wald test was carried out to determine the statistical difference between the genetic parameter estimates and zero. In addition, correlations between the genomic breeding values for GESS of the reproduction traits, as well as production, health, and type traits, were estimated using the method proposed by Calo et al. (1973) and bull EBVs (filtered based on EBV reliability, as described below). Standard errors (SE) of the approximate genetic correlations were calculated based on the formula proposed by Sokal and Rohlf (1981). The predicted transmitting ability (PTA) of six production traits (milk yield, milk protein yield, milk protein ratio, milk fat yield, milk fat ratio, and somatic cell score), two health traits (reproductive disease and udder disease), and five type traits (dairy character, milking system, capacity, rump, and overall conformation) with reliabilities higher than 20, 20, and 30%, respectively, were provided by BDCC. The statistical models used for the genetic evaluation of these traits are described in Miglior et al. (2009) and Wu et al. (2013). Individuals with PTA reliabilities for GESS on reproductive traits above 40% were used for the calculation of correlation of breeding values.

Genome-Wide Association Study and Functional Enrichment Analyses
The "Fixed and random model Circulating Probability Unification" (FarmCPU) R package  was used to perform single-SNP regression analyses. FarmCPU is a multi-locus model that incorporates multiple markers simultaneously as covariates to partially remove the confounding effect between testing markers and kinship . De-regressed proofs (DRP) of the GESS and GED for female reproductive traits were derived following the procedures suggested by Wiggans et al. (2011). Individuals with accuracies greater than 10% were used as dependent variables in the GWAS model. The obtained p-values were corrected for multiple testing by calculating the false discovery rate (FDR) (Benjamini et al., 2001) at the 5% genome-wise level. Quantile-quantile (Q-Q) plots and the inflation factor λ (Devlin and Roeder, 1999) were used to investigate population stratification by comparing the observed and expected distributions of -log(p-value).
Positional genes located at up to 200 kb upstream and downstream of the significant SNPs were identified using the BiomaRt package (Durinck et al., 2005). This 200-kb window was defined based on the linkage disequilibrium level of the studied population (Supplementary Figure 1). The ClueGO module in Cytoscape (Bindea et al., 2009) was used to identify candidate genes, Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms. Furthermore, the Cattle Quantitative Trait Locus (QTL) database (Cattle QTLdb Release 42 3 ) was used to identify important trait-QTL associations previously reported in the literature.

Descriptive Statistics of Phenotype
The descriptive statistics of the five reproductive traits evaluated are shown in Table 1. For the insemination success traits (CR and NRR56), the mean CR was 43 ± 49% and the mean NRR56 was 50 ± 50% in Chinese Holstein cattle. These estimates are lower than those reported by Guo et al. (2014), which were inferred from the number of services recorded from 2001 to 2010. The average CE, SB, and GL were 1.06 ± 0.27, 1.07 ± 0.25, and 278.36 ± 6.18 days, respectively. The proportion of CE scores higher than 1 was 5.21%, and the SB rate was 6.67%. Vieira-Neto et al. (2017) reported a mean GL of 276 ± 6 days in Holstein cows, which is slightly lower than that of the current study. GL is affected by various environmental effects, such as temperature (McClintock et al., 2003) and cow parity number (King et al., 1985), which may have led to the discrepancy of GL between different populations.

Heritability and Repeatability Estimates
The (co)variance components, heritability, and genetic correlations for GESS and GED for all five reproduction traits are shown in Table 2. Overall, the heritability estimates for the success traits were low, but the heritability of GESS were significantly higher than those of GED (0.020 ± 0.004 vs. 0.011 ± 0.000 for CR and 0.024 ± 0.005 vs. 0.013 ± 0.001 for NRR56), indicating that service sires actually have greater genetic impact on insemination success than the mating cows. Fertility traits are known to have low heritability (VanRaden et al., 2004), which makes it more difficult to obtain faster genetic progress compared to more heritable traits. Therefore, larger datasets and novel phenotypes (Fleming et al., 2019) should be generated for increasing genetic progress for fertility performance. The small, but significant, genetic contribution of service sires on success traits indicates the possibility of genetically improving the GESS. The repeatability estimates of these reproduction traits were almost equal to their heritabilities, except for the GED on NRR56 (0.031 ± 0.001 vs. 0.013 ± 0.001), indicating the inconsistency among these records. Hence, more repeated records are needed for greater EBV accuracies. Another reason for the low repeatability might be the ignorance of non-additive genetic effects, which could lead to the underestimation of genetic parameters. For most traits, the genetic variance of the GESS was comparable to that of the GED, demonstrating that service sires have considerable genetic contribution to female reproductive outcomes. The heritability estimates for female CR in other studies range from 0.01 to 0.03 (Azzam et al., 1988;Bagnato and Oltenacu, 1993;Muuttoranta et al., 2019). These low estimates support the estimates of GED on CR in this study. The heritability of GED on NRR56 was consistent with the value (0.01) reported by Sun and Su (2010) but different from that of Hoekstra et al. (1994) (0.04) and Eghbalsaied (2011) (0.002-0.003). These discrepancies may be attributed to the differences in the populations evaluated and the statistical models used since GESS was only included in the current study. Some studies have considered the GESS as a non-genetic random effect apart from the GED, and the heritabilities obtained from such models were shown to be lower than 0.02 for CR and NRR56 (Weigel and Rekaya, 2000;Kuhn and Hutchison, 2008). However, Tiezzi et al. (2011) suggested that the heritabilities generated using the approach mentioned above are lower than those fitting GESS. Furthermore, Tiezzi et al. (2013) reported that the heritability and repeatability of the GESS on NRR56 were both 0.01 for Italian Brown Swiss cattle. The characteristics of the two populations evaluated might have caused the discrepancy in the estimates observed. In the study of Tiezzi et al. (2013), error of heritability of dam; re d , repeatability of dam; SE re d , standard error of repeatability of dam; r, genetic correlation between genetic effect of service sire (GESS) and genetic effects of the dam (GED); SE r , standard error of genetic correlation between GESS and GED; CR, conception rate; NRR56, 56-day non-return rate; CE, calving ease; SB, stillbirth; GL, gestation length.
CR, conception rate; NRR56, 56-day non-return rate; CE, calving ease; SB, stillbirth; GL, gestation length. a The EBVs of 400-500 individuals with reliability greater than 40% for reproductive traits and 20% for other traits were used for calculating the approximate correlations. b The EBV of milk fat ratio and milk protein ratio, respectively, come from the formulas below: EBV fat ratio = EBV fat yield × 100 − EBV milk yield × fat ratio EBV milk yield + milk yield EBV protein ratio = EBV protein yield × 100 − EBV milk yield × protein ratio EBV milk yield + milk yield where EBV fat ratio , EBV fat yield , EBV milk yield , EBV protein ratio , and EBV protein yield represent the EBV of milk fat ratio, milk fat yield, milk yield, milk protein ratio, and milk protein yield, respectively; fat ratio , protein ratio , and milk yield represent the average milk fat ratio, milk protein ratio, and milk yield of cows in their second lactation.
Frontiers in Genetics | www.frontiersin.org Frontiers in Genetics | www.frontiersin.org a smaller and older population than that in the current study was used, and the cattle in their population performed better in NRR56 (0.70 for average), indicating a better management and genetic level of their herds.
For calving traits, the heritability estimates of GESS are lower than the GED for CE (0.004 ± 0.001 vs. 0.029 ± 0.002) and similar for GL (0.102 ± 0.001 vs. 0.102 ± 0.000) but higher for SB (0.078 ± 0.010 vs. 0.008 ± 0.001). Both GESS and GED show   low heritability for SB and CE, though a moderate estimation was observed for GL. Except for the GED for SB (0.021 ± 0.002) and GL (0.116 ± 0.004), the repeatabilities were almost equal to the corresponding heritabilities. Each cow has only one calving record in each parity, which results in poor consistency among repeated data. The low heritability estimates of SB and CE indicate that these traits might benefit more from genomic information. A S-MGS model was used to evaluate calving traits in previous studies, which resulted in low heritability for CE and SB (0.07-0.13) Heringstad et al., 2007). The S-MGS model converted the direct and maternal variances to sire and maternal grandsire variances, which was, to some extent, different from those of the present study. Due to the differences in the model used in the current study, here we simply compare our findings with the results of the S-MGS model from previous studies. The maternal heritabilities estimated by other models ranged from 0.01 to 0.13 for GL (Jamrozik et al., 2005;Crews, 2006;Mujibi and Crews, 2009) and 0.02-0.08 for CE and SB (Luo et al., 2002;Jamrozik et al., 2005;Eaglen et al., 2013). The GED variance obtained in the current study for SB was 0.0002, which is much lower than the estimates from a previous study with Norwegian Red cows (Heringstad et al., 2007). The heritability of the GED was higher than that of the GESS for CE, which may have been caused by a larger maternal effect (Jamrozik et al., 2005)-for instance, cow body conformation is genetically associated with calving difficulty (Dadati et al., 1985;Ga et al., 2021), indicating the crucial role of GED on CE.

Genetic Correlations
Former genomic analyses reported that some genes (e.g., SPP1) regulate both tissue and embryonic growth (Weintraub et al., 2004;Rangaswami et al., 2006), which is important for both the male and female aspects of calving traits. Furthermore, some genes related to spermatogenesis have been shown to affect cow reproduction performance (Peddinti et al., 2010;Li et al., 2012b;Buzanskas et al., 2017). Thus, the genetic covariances of male and female were considered. The genetic correlations between the GESS and GED were significantly different from zero for SB and GL (P < 0.05; based on a t-test). The correlation coefficient was considerably high in GL, whereas that in CE was low and non-significant (P > 0.05; CE: 0.099 ± 0.081; SB: 0.139 ± 0.071; GL: 0.787 ± 0.040). The results indicate the homogeneity between the GESS and the GED for SB and GL. The genetic covariances of the insemination traits between the effects of service bull and dams were usually ignored in previous studies (Berry et al., 2011). We evaluated these covariances because there are evidence of low but statistically significant correlations between GESS and GED [e.g., NRR56: 0.010 ± 0.002; (Tiezzi et al., 2013)]. However, there was no significant correlation between the two terms in the current study (CR: 0.044 ± 0.066 and NRR56: −0.094 ± 0.063). The genetic correlations between the direct and maternal effects for SB and CE can be quite variable, ranging from −0.24 to 0.12 (Luo et al., 1999;Steinbock et al., 2003;Wiggans et al., 2003;Cole et al., 2007;Heringstad et al., 2007;Vanderick et al., 2014). For GL, the correlations are usually negative and stronger (−0.13 ∼ −0.85) (Cubas et al., 1991;Bennett and Gregory, 2001;Hansen et al., 2004;Crews, 2006;Cervantes et al., 2010). These discrepancies are reasonable, because we considered both the direct effect of dam and maternal effect as dam effect. The correlations of breeding values of the GESS of reproductive traits as well as production, health, and type traits are shown in Table 3. For all the traits, the genetic information of approximately 400-500 individuals with reliability greater than 40% for reproductive traits and 20% for the other traits was used for the calculation of breeding value correlations. Interestingly, the GESS was negatively correlated with most production and type traits (e.g., milk yield and overall conformation), while positive correlations were observed between GESS and health traits such as reproductive disorders. We found that the GESS on CR was unfavorably and significantly related to milk yield (−0.218 ± 0.039), indicating that selection exclusively on milk production might indirectly result in a decline of insemination performance of the service sire. Murray et al. (1977) reported a negative correlation between male fertility and milk yield FIGURE 2 | Manhattan plots of the genome-wide association studies for genetic effect of service sire on reproductive traits. The x-axis and the y-axis represent the chromosome number and the observed −log 10 (P−value), respectively. The single-nucleotide polymorphisms were plotted against their genomic positions. The lines in the plots indicate the thresholds of false discovery rate (0.05) in the corresponding traits: (A) conception rate, (B) 56-day non-return rate, (C) calving ease, (D) stillbirth, and (E) gestation length.
(−0.26), which is in contrast to the positive correlation (0.13-0.29) reported by Raheja et al. (1989). Further verification about the biological correlation between male fertility and milk production is needed due to the inconsistent relationships observed. The genetic correlations between GESS on success traits and reproductive disease were positive (0.174 ± 0.043 for CR and 0.203 ± 0.042 for NRR56). Progesterone is regarded as a responsible factor for cattle ovarian follicular cysts (Silvia et al., 2002). Ramal-Sanchez et al. (2020) also suggested that progesterone is related to sperm release, which affects the fertilization ability of spermatozoa. Therefore, reproductionrelated hormones might account for the potential biological relation between GESS on success traits and reproductive disease. Cows with a higher incidence of ovarian cysts tend to have lower fertility. The genetic correlation between GESS on GL and overall conformation was negative, indicating that undesirable body conformation might lead to longer GL, with worse development after birth (Bourdon and Brinks, 1982). These correlations indicate that direct selection for production, health, and type traits may have a favorable effect on service sire calving performance but may lead to an unfavorable decline in service sire mating performance.

Genome-Wide Association Studies
The estimated breeding values obtained from former genetic estimation with accuracy above 10% were used for DRP calculation. Hence, the EBV of 2,996 and 1,147 individuals were used in GWAS for GESS and GED, respectively. The detailed information of the genomic regions and candidate genes for GESS and GED for the five reproductive traits are summarized in Tables 4, 5. In addition, the Q-Q plots and Manhattan plots are provided in Figures 1-4. The Q-Q plots and λ of GESS on 56-day non-return rate indicated a slight inflation of the results. A total of 162 significant SNPs were detected, including two significant SNPs located in the X chromosome (pseudoautosomal region) (Johnson et al., 2019). The P-values ranged from 1.34 × 10 −10 to 1.05 × 10 −5 , and FDR ranged from <0.001 to 0.050. Furthermore, we mapped the significant SNPs to the Cattle QTL database (see Text Footnote 3), and the overlapped QTL regions are listed in Tables 6, 7. For GESS, 100 SNPs were found to be significant, with 23, 24, 18, 18, and 17 associated with CR, NRR56, CE, SB, and GL, respectively ( Table 4). For the success traits, 53 nearby annotated protein-coding genes located 200 kb upstream and downstream of the significant SNPs were mapped, including genes related to sperm development (e.g., BMP4) and early embryogenesis (e.g., LRRC34). Han and Peñagaricano (2016) also identified a genomic region (1.5 Mb) located on BTA13 which explained more than 0.50% of the total additive genetic variance for SCR (male fertility). This region contains the CTCFL gene detected in the present study. Although previous studies have attempted to identify candidate genes related to GESS on success traits (Taylor et al., 2018), novel candidate genes were identified in this study-for instance, the gene BMP4 (bone morphogenetic protein 4) was previously indicated as a candidate gene for spermatogenesis (Hu et al., 2004;Li et al., 2014) and follicle development (Nilsson and Skinner, 2003;Shimizu et al., 2004;Fatehi et al., 2005;Li and Ge, 2011). These associations indicate the paternal and maternal effects on the pre-implantation stage  of embryo development, respectively (Koide et al., 2009;Li et al., 2012a). In this context, the present study also revealed the GESS on female conception. Ou et al. (2020) carried out a transcriptome analysis on spermatogonial stem cells (SSCs) and reported that LRRC34 was highly expressed in mouse SSCs and was essential for in vitro SSC proliferation. RXFP3 was identified to be differentially expressed in human sperm and was likely diminished in spermiogenesis (Heidari et al., 2018). Furthermore, other studies suggested PPP2R2B and PCK1 as candidate genes affecting the semen quality traits in livestock Gao et al., 2019).
One hundred fourteen nearby proteinase genes located 200 kb upstream and downstream of the significant SNPs related to the GESS of calving traits were mapped. The fetal growth and metabolism show a specific pattern during pregnancy (Bell et al., 1993) and supposedly reflect on calving traits. Some potential genes related to calving traits in our study were previously associated with carcass and meat quality traits, including MTUS1 (Albrecht et al., 2016), PLCH1 (Lemos et al., 2016), F2RL1 (Srikanth et al., 2020;Zhang et al., 2017), MYO7B (Doyle et al., 2020;Jia et al., 2020), WWOX (Grigoletto et al., 2020), TFB2M (Jiang et al., 2006;Song et al., 2019), andSMYD3 (De Vos, 2018). Furthermore, F2RL1 was reported to affect the body size in Chinese Holstein cattle (Zhang et al., 2017). Some other genes were identified as essential genes for embryogenesis [SEMA4D (Masuda et al., 2004), and CCNG2 (Ma et al., 2015) and CKS2 (Martinsson-Ahlzén et al., 2008), which contribute to subsequent fetal development]. Perkins et al. (1995) collected human fetal plasma samples both at mid-gestation and parturient to explore the trend of corticotrophin-releasing hormone-binding protein (CRHBP) during the different gestation stages and reported CRHBP as being functional in both maternal and fetal circulation. F2RL2, LIMS2, and LIMS1 were indicated by Forde et al. (2012) as pregnancy-associated genes that are differently expressed in the endometrium of cattle during early pregnancy, indicating the potential role of these genes on successful pregnancy establishment. The SNP rs43354413 located in BTA 3 was identified as significant for GESS on NRR56 and GL, with 10 protein-coding genes in close proximity-for instance, the SSBP3 gene (located approximately 131 kb upstream of this marker) was reported to regulate mouse embryonic stem cell differentiation to trophoblast-like cells (Liu J. et al., 2016). These findings suggest the possible function of SSPB3 of GESS on reproduction performance. We used Chinese Holstein population and redefined the genetic components of calving traits in current study, but some SNPs (rs42813960 and rs110402487) in our study are consistent with previous genomic studies of calving traits focused on BTA18 (Müller et al., 2017), indicating the potential importance of this chromosome in calving performance. Fang et al. (2019) proposed ZNF613 as a candidate gene for paternal contributions to GL, and in our population, we detected one GLrelated marker located downstream of ZNF613 (rs110402487). Furthermore, we also mapped these SNPs to the Animal QTL Database (see Text Footnote 3) and subsequently found that some of them were located in QTLs associated with related traits (Table 6). Particularly, rs108993952 was related to GESS on GL and meanwhile located in the QTL related to some maternal calving traits, which is a promising candidate marker for calving performance.
Sixty-two SNPs were found to be significant for GED,with 13,12,9,17,and 11 being associated with CR, NRR56, CE, SB, and GL, respectively (Table 5). One hundred sixtyseven protein-coding genes were found within 200 kb of those SNPs. rs136577145 was significant for both GESS and GED on GL and was about 63 kb upstream of rs41578821, which was detected as significant for GED on CE. The nearby genes mapped with rs136577145 were TFB2M, CNST, SCCPDH, H3-5, and SMYD3. TFB2M, a mitochondrial transcription specificity factor, as well as SMYD3, a histone lysine methyltransferase, were previously linked to bone and skeletal muscle tissue development (Norrbom et al., 2010;Fujii et al., 2011;Sun et al., 2015)-that is, both paternaland maternal-derived effects on GL were coincident with growth ability. In addition, WWOX overlapped between GESS and GED on other calving traits, though no gene overlapped between GESS and GED on success traits (Figure 5). This is consistent with the genetic correlation estimates observed ( Table 2). The low overlap between GESS and GED on success traits demonstrates that different genes are involved in the insemination outcomes from the dam and service sire contributions.   Except for the genes mentioned before, the genomic results of GED are consistent with previous studies (e.g., SMAD4). For GED on success traits, 86 nearby protein-coding genes located 200 kb upstream and downstream of significant SNPs were mapped. Particularly, SMAD4 has been reported to be linked to embryogenesis and folliculogenesis in mice (Chu et al., 2004;Pangas et al., 2006). Lee et al. (2014) reported SMAD4 to be involved in early embryonic development in cattle and to regulate the effects of follistatin, which influences the environment-independent part of the interval from the first to the last insemination in cattle . Eighty-one protein-coding genes were mapped for GED on calving traits and previously considered as essential genes for the maternal mechanism during pregnancy, including CXCL10 (Walker et al., 2012), HBEGF (Jessmon et al., 2009), PCDHA13 (Lotfan et al., 2018), TGM2 (Purfield et al., 2019), C1QB (Cochran et al., 2013), CYSTM1 (Purfield et al., 2019), and EPHB2 (Purfield et al., 2015), which were also identified as candidate genes for maternal effect on calving traits in dairy and beef cattle. We also observed that some significant SNPs located in QTLs are associated with related traits (Table 7)-for instance, rs42427669, which was related to GED on NRR56, is located in a QTL region associated with CR (Kiser et al., 2019), while the GL-related marker rs108993952 is located in QTLs associated with CE and SB .

Functional Enrichment Analyses
The enriched GO terms and KEGG pathways that passed the criteria of the Benjamini-Hochberg-corrected p-value < 0.05 are summarized in the Supplementary Files, and the genes shared between the main terms or pathways are presented in Supplementary Figures 2-5. Noticeably, both GESS and GED on success traits were enriched in neural development-related terms, such as neural crest cell migration (GO:0001755) for GESS on success traits and positive regulation of neuron projection development (GO:0010976) for GED on success traits. For GESS on calving traits, genes were enriched mainly in the thrombinactivated receptor signaling pathway (GO:0070493), microvillus (GO:0005902), neural precursor cell proliferation (GO:0061351), liver development (GO:0001889), N-methyltransferase activity (GO:0008170), sulfotransferase activity (GO:0008146), and cyclin-dependent protein kinase holoenzyme complex (GO:0000307). These categories include functionable genes named CCNG2 and CKS2 (GO:0000307), LIMS2 (GO:0001889 and GO:0061351), MYO7B, WWOX (GO:0005902), SMYD3, TFB2M (GO:0008170), F2RL1, and F2RL2 (GO:0070493) as previously discussed. In addition, the enrichment analysis for GED showed main terms such as CXCR3 chemokine receptor binding (GO:0048248), homophilic cell adhesion via plasma membrane adhesion molecules (GO:0007156), and synapse pruning (GO:0098883), including CXCL10 (GO:0048248), PCDHA13 (GO:0007156), and C1QB (GO:0098883). Therefore, our findings provide further evidence of the possible genetic mechanism of GESS and GED on reproductive performance in Holstein cattle.

Future Prospects
Generally, only the GED for reproduction performance is analyzed in genetic evaluations, though research including the current study have shown the need to dissect GESS (Tiezzi et al., 2011). Compared to previous studies, we fitted an improved model considering the covariance between GESS and GED and also identified candidate genes associated with GESS and GED. The GESS on reproductive traits is small but significant. The low repeatability estimates indicate the poor consistency among repeated records. Therefore, more accurate records and novel traits are required. With recent improvements in data collection, GESS might become an important factor on the genetic evaluation of reproductive performance. Genomic selection is also expected to contribute to improve the accuracy of breeding value for these lowly heritable traits (Rice and Lipka, 2019).
Additional analyses with larger datasets and in independent populations (e.g., different breeds) are recommended.

CONCLUSION
The GESS on reproductive traits is heritable, with a similar genetic variance to the GED. Moreover, the approximate genetic correlation among the GESS and production, health, and type traits is unfavorable for the success traits (CR and NRR56) but favorable for the calving traits (CE, SB, and GL). A total of 100 and 62 significant SNPs were detected to be associated with GESS and GED on those five reproductive traits, respectively. Among them, five genes (BMP4 and CTCFL for success traits and WWOX, TFB2M, and SMYD3 for calving traits) are suggested as important candidate genes for GESS according to positional and functional analyses. As GESS and GED are lowly heritable, genomic prediction might be a promising alternative for breeding schemes aiming to improve fertility performance in dairy cattle.

DATA AVAILABILITY STATEMENT
The data analyzed in this study is subject to the following licenses/restrictions: this manuscript utilizes proprietary data.
Requests to access these datasets should be directed to YW, wangyachun@cau.edu.cn.

ETHICS STATEMENT
The studies involving animals were reviewed and approved by the Animal Welfare Committee of the China Agricultural University (Protocol Number: DK996).

AUTHOR CONTRIBUTIONS
ZC and YW conceived the study. ZC and LB designed the analyses, discussed the results, and drafted the manuscript. ZC performed all the data analyses and the data curation of genotype. ZC, HL, and RS prepared the phenotypic files. YC assisted with the modeling analyses. LL and GG provided support for the collection of raw data. All authors have read and approved the final version of the manuscript.

FUNDING
This study was supported by the China Agriculture Research System of MOF and MARA, the program for Changjiang Scholar and Innovation Research Team in University (IRT_15R62), and the National Agricultural Genetic Improvement Program (2130135). This study also received funding from Beijing Sanyuan Breeding Technology Ltd., Co. funded project (SYZYZ20190005). The funder had the following involvement with the study: study design. All authors declare no other competing interests.

ACKNOWLEDGMENTS
The authors acknowledge Beijing Sanyuan Breeding Technology Company Limited for providing access to the datasets used. The authors also grateful to Guosheng Su from Aarhus University for providing suggestions about the statistical model of the current study and to Qing Xu from Beijing Jiaotong University for giving advice on the genomic analyses.