Health and breeding records for 5,750 primiparous Holstein cows were obtained using Dairy Comp 305 (Valley Agricultural Software, Tulare, CA, USA) records from a single dairy in southern Idaho. Data were recorded daily for breeding, production, and disease. The dairy consists of a milking herd of 2,300 cows housed in a dry lot with shade and fed a total mixed ration. Only cows that were pregnant to the first AI (n = 748) or were confirmed pregnant to the first ET (n = 277) were included in the study to eliminate possible confounding of phenotypes for embryonic loss. Cows were bred following a double ovulation synchronization protocol or after standing estrus prior to AI. There was no difference (p < 0.05, ANOVA) in the frequency of SA in cows bred by AI during observed estrus or following synchronization. Cows that were ET recipients received a double ovulation synchronization protocol with embryo transfer on day 7 post-estrus. Pregnancy was determined using trans-rectal ultrasound at 30 days post-breeding. Phenotypes were based on whether a cow maintained a pregnancy to ≥ 260 days of gestation (control) or spontaneously aborted between days 42 and 260 of gestation (case). Cases were identified by the cows returning to estrus following a confirmed pregnancy. AI controls were compared to AI cases, while ET controls were compared to ET cases. Four AI cows and three ET recipient cows were removed from the study as they experienced mastitis, metritis, metabolic issues, lameness, or respiratory disease between approximately day 30 before breeding and/or before calving. After filtering for these events, 679 AI controls, 69 AI cases, 236 ET controls and 41 ET cases remained in the study. All cows were at least 25 months of age at the beginning of their first lactation.

Genotypes from Zoetis (Kalamazoo, MI, USA) were imputed to a higher density utilizing Beagle (V. 4.1) and an in-house reference population of 4,800 US Holsteins genotyped with the Illumina (San Diego, CA, USA) BovineHD BeadChip (15). The 636,042 imputed single nucleotide polymorphisms (SNPs) had a call accuracy of >95% across the genome. Accuracy was determined by reducing the genotyping density of animals in the reference population to approximately 50,000 SNPs. Imputation was then performed on the lower density genotypes to restore them to the BovineHD BeadChip level for assessing the accuracy of imputation. Genotypes from the imputed SNPs were compared with the SNPs obtained from the initial genotyping using the BovineHD BeadChip.

Prior to the GWAA, imputed genotypes were filtered for quality. SNPs were removed if they had a call rate <0.90, had a minor allele frequency (MAF) <0.01 or if they failed (p < 1 x 10⁻¹⁰⁰) Hardy–Weinberg equilibrium testing. Primiparous cows bred by AI had 4,733 SNPs removed for call rate, 140,423 SNPs were removed for MAF, and 1,012 SNPs were removed for failing Hardy–Weinberg equilibrium testing. For primiparous cows receiving ET 4,658 SNPs were removed for call rate, and 151,385 SNPs were removed for low MAF. No SNPs failed Hardy–Weinberg equilibrium testing for cows that were ET recipients, nor were any cows removed for genotyping call rate (<0.90) or for differences between genotypic and phenotypic sex. After quality control was complete, a total of 1,025 cows remained for the analyses, with 489,927 SNPs remaining for cows bred by AI and 479,927 SNPs remaining for ET recipient cows.

A GWAA was performed using the SNP and Variation Suite (SVS) software version 8.1 (70). An Efficient Mixed Model Associated eXpedited (EMMAX) statistical approach was used with an identity by state matrix (16). The general mixed model for EMMAX is expressed as γ = Xβ+Zμ+ε where y = the vector of observed phenotypic values, X = a matrix of fixed effects, β = regression coefficients, Z = a matrix containing the observed random effects, u = vector of random effects related to variants of allele substitutions in the population, and ε = residual effects (16). Associations with SA were identified when FDR <0.05 (17, 18). Dosage compensation for SNPs on the X chromosome were not accounted for in EMMAX.

Due to the unknown inheritance of SA in cattle, three inheritance models (additive, dominant, and recessive) were analyzed for AI and ET groups. The additive model identified associations where having two minor alleles (aa) was twice as likely to impact the phenotype (SA) compared to having no minor alleles (AA), and half as likely to have an impact when one minor allele (Aa) was present. The dominant model identified associations where having one or two minor alleles (Aa or aa) was likely to have more impact on the phenotype compared to not having a minor allele (AA). A recessive model identified associations where having two minor alleles (aa) was likely to have more impact on the phenotype compared to having one or no minor allele (Aa or AA).

A principal component analysis was performed to evaluate population stratification. Two distinct clusters within the population based on birth year were identified (Supplementary Figure S1). An ANOVA (p < 0.05) was used to test if birth year, service sire, and season of conception affected SA. The tested covariates with the associated p-values for each population are listed in Table 1. To ensure that the population stratification was accounted for, the genomic inflation factor lambda (λ_GC) was calculated in SVS for each model. For AI-bred primiparous cows, no covariates were included as none of the tested covariates showed a significant (p < 0.05) effect on SA (Number of Sires used). In contrast, ET recipient primiparous cows had a difference (p < 0.05) in SA by season, therefore, this covariate was included in the analysis. Additional covariates, not identified by ANOVA as significant, were not included in the analysis to minimize false negative errors due to over-correction of possible population stratification.

If multiple SNPs were associated with SA on a chromosome, loci were defined by a threshold of D'> 0.7 (19, 20). For each locus, positional candidate genes were identified within the average haplotype size of 30.5 kb, calculated using the method of Gabriel (21), 5′ or 3′ to the associated SNP based on the bovine ARS-UCD 1.2 genome assembly (22).

The proportion of variance explained by a SNP was calculated in SVS utilizing the notation of Further Optimization when Covariates are Present, as described by Segura (72) and Vilhjalmsson (73). The sum of the proportion of variance explained for all SNPs will exceed 100% due to SNPs within a locus not being independent. This calculation was used to evaluate the percentage of total variation observed in spontaneous abortion that can be attributed to these SNPs.

A measure of the strength of the association between having the risk allele and experiencing a SA, referred to as relative risk, was calculated utilizing the frequencies of the minor and major alleles (23). The following equations were utilized for calculating relative risk for SA: a(a+b)/c(c+d), where a is the minor allele frequency in cases, b is the minor allele frequency in controls, c is the major allele frequency in cases, and d is the major allele frequency in controls.

Heritability was estimated for SA using a genomic best linear unbiased predictor (GBLUP) analysis (24, 25) with the average information algorithm (AI-REML), which is a bivariate restricted maximum likelihood analysis (26, 71). The AI-REML GBLUP method calculates variance components that are then used to calculate heritability, as the pseudo-heritability estimated by EMMAX can over-inflate the estimate in limited sample sizes.

Literature and database searches were conducted to identify any shared genomic regions or genes associated with SA and other fertility and production traits. These traits included the fertility index in Nordic Red Cattle, which encompasses multiple traits such as the number of inseminations per conception, duration of the interval from calving to first insemination, days from first to last insemination, and 46-day non-return rate (27). Other fertility traits included cow conception rate (percentage of cows that are pregnant at each service), daughter pregnancy rate (percentage of cows that become pregnant during each 21-day period), days open, and early embryonic loss (20, 28–31). Additionally, a comparison was made using tissue expression data based on fertility status and expression in reproductive tissues (32, 33). Any shared loci or positional candidate genes will be listed and discussed in the Discussion section.

Additionally, the comparison of production traits was conducted in two ways: first, through a literature search for shared loci and positional candidate genes related to production traits such as milk yield and milk components like fat yield, fat percentage, protein yield, and protein percentage (34). Secondly, multiple GWAA were utilized on the population of animals in this study that had complete lactation records for the following three production traits: total milk, total fat, and total protein.

The highest incidence of SA for primiparous cows bred by AI occurred between days 42 and 91 of gestation. The single highest day for SA was day 45 of gestation, with 12 cows aborting (Figure 1). These gestational days when SA was noted represent the day the cow returned to estrus after being identified as pregnant by trans-rectal ultrasound at approximately day 30 of gestation. The 54 SA occurrences between days 42 and 91 of gestation accounted for 78% of the fetal loss in this population. The second greatest interval of gestation associated with SA was between days 168 and 260, during which 20% of fetal loss occurred.

The λ_GC for SA in AI cows was 0.99 for the additive inheritance model, 1.01 for the dominant inheritance model, and 0.96 for the recessive model. The estimated heritability for SA in primiparous cows bred by AI was 0.03 (±0.053). As the standard error for the heritability estimate is larger than the estimate itself, this provides low confidence in the estimated heritability. Since the majority of the loci associated with SA were recessive, and heritability measures additive genetic variance, it is unsurprising that the heritability estimate for SA is low. This estimated heritability for SA is similar to that reported for primiparous cows at 0.09 (±0.08) in a previous study (3) and for other reproductive traits such as the number of services, days from calving to first service, days open, and calving interval, which have been reported to be 0.040 ± 0.017, 0.034 ± 0.011, 0.053 ± 0.019, and 0.056 ± 0.014, respectively (3, 37). Rahbar et al. (38) also reported similar heritability estimates for several fertility traits, such as days open (0.016) and gestation length (0.123). Two other traits with reported heritability estimates comparable to that calculated in this study include the interval from first calving to first insemination (0.013) and the interval from first to successful insemination for primiparous cows (0.038) (39).

In the recessive model, there were 85 loci associated (FDR <0.05) with SA (Figure 2, Supplementary Table S1). The most significant loci associated with SA in the recessive model were on BTA3 at 110 Mb (FDR = 2.9 × 10⁻⁵), BTA22 at 3 Mb (FDR = 4.7 × 10⁻⁵), BTA24 at 37 Mb (FDR = 3.5 × 10⁻⁵), and BTA28 at 9 Mb (FDR = 1 × 10⁻⁴) and 39 Mb (FDR = 3.07 × 10⁻⁵). Positional candidate genes at these loci include: SFPQ, ZMYM4, KIAA0319L, RBMS3, DLGAP1, MTR, TRNAE-UUC, LOC112444209, LOC112446145, LOC618220, and LOC107131937. An additional 156 positional candidate genes were identified as associated with SA, for a total of 168 positional candidate genes associated with SA in primiparous cows bred by AI. The positional candidate genes were grouped into 16 functional groups as shown in Supplementary Table S2.

The highest proportion of phenotypic variance explained for the primiparous cows bred by AI was 5.6% from the locus on BTA28 at 9 Mb, with two positional candidate genes, MTR and TRNAE-UUC. Additionally, 21 loci representing 45 positional candidate genes, listed in Table 2, had a proportion of variance explained >3%. The average relative risk for the most significant SNPs associated with SA in primiparous cows bred by AI was 2.20 (Supplementary Table S3). The four positional candidate genes of the most significant loci are further described below.

Positional candidate genes associated with SA in primiparous cows bred by AI were differentially expressed in the endometrium and in different cell types of the bovine placenta (Supplementary Figure S2) during the embryonic and fetal periods. Twenty-eight positional candidate genes were differentially expressed in high fertility compared to sub-fertile heifers at day 17 of gestation (32). An additional 42 positional candidate genes were expressed in the reproductive tissues of cattle, based on single-cell data during the embryonic and fetal periods, which support their involvement in pregnancy (Table 3; Supplementary Figure S2) (33).

Primiparous cows that were ET recipients also experienced the greatest number of SA on day 45 of gestation (Figure 3). Seventy-six percent of the SA in ET recipients occurred during the first trimester.

The λ_GC for SA in cows that were ET recipients was 1.02, 1.02, and 1.03 for the additive, dominant, and recessive inheritance models, respectively. No loci were associated with SA in cows that were ET recipients in the additive or dominant inheritance models, but four loci were associated (FDR <0.05) with SA in the recessive inheritance model on BTA9 and BTA13 (Figure 4, Table 4). There were 14 positional candidate genes within the two loci on BTA13 (Table 4). All four loci associated with SA in ET recipient primiparous cows had a proportion of variance explained >3%. The smallest variance explained was 8.7%, belonging to the locus on BT 13 at 46 Mb, while the largest proportion of variance explained was 10.5% for the locus on BTA13 at 47 Mb. The positional candidate genes for these loci can be seen listed in Table 4.

Loci and positional candidate genes associated with SA were also compared with production traits to identify possible correlations to determine if selection for improved fertility would negatively impact selection for production traits. The loci and positional candidate genes associated with SA in AI bred and ET recipient primiparous cows were not associated (p <0.05) with dairy production traits. These traits included total milk, total fat, total protein, milk yield, fat yield, fat percentage, protein yield and protein percentage as mentioned in the materials and methods. This indicates that selection for production traits while selecting for reduced SA could be done simultaneously without reducing genetic progress in production traits.

One locus on BTA17 at 51 Mb was associated (FDR = 0.049) with SA in the additive model with the positional candidate gene nuclear receptor corepressor 2 (NCOR2) (Figure 2, Supplementary Table S1). The NCOR2 gene codes for a nuclear corepressor that interacts with multiple transcription factors. This corepressor orchestrates proliferation and differentiation during B cell development and is critical for fetal development (44). There were six loci associated with SA in cows bred by AI that were previously identified as associated with fertility traits. For a locus to be determined to be shared, the associated SNPs were required to be within a haplotype (31 kb) of each other. Two of these six loci, located on BTA1 and BTA4, were shared with genomic regions identified by Cole et al. (28) for daughter pregnancy rate (DPR) in Holstein cows. The sharing of loci between DPR and SA may be due to similar physiological processes necessary for the initiation and maintenance of pregnancy in cattle. A locus on BTA16 associated with cow conception rate to the first service and number of times bred in Holstein cattle was also associated with SA (31). A study by Höglund et al. (45) identified two loci (BTA2 at 46 Mb and BTA29 at 48 Mb) that were associated with fertility traits in Nordic cattle and shared with SA in this study. The locus on BTA2 was also shared with a genomic region identified for the Nordic female fertility index and days from first to last insemination in Nordic Holstein, Nordic Red, and Jersey cows (45). The Nordic female fertility index is a multi-trait index based on data collected on dairy cattle in Denmark, Sweden, and Finland and includes the number of inseminations required per pregnancy in cows and heifers, the length in days of the interval from calving to first insemination in cows, days from first to last insemination in cows and heifers, 56-day non-return rate in cows and heifers, and heat strength in cows and heifers (45). A locus on BTA12 at 52 Mb was also associated with a fertility index in Nordic Red cattle and with SA in Holstein cattle (27). Finally, a locus on BTA29 that was associated with SA was also associated with the number of inseminations per conception and days from first to last insemination in cows and heifers in Nordic Holstein, Nordic Red, and Jersey cows (45). Additionally, there were nine loci (BTA1 at 128 Mb, BTA5 at 85 Mb, BTA7 at 6 Mb, BTA7 at 85 Mb, two on BTA13 at 75 Mb, BTA17 at 35 Mb, BTA18 at 9 Mb, and BTA23 at 43 Mb) associated with SA in heifers bred by AI (13) and primiparous cows bred by AI. This indicates that shared genomic regions across parities could be used for genomic selection to help lower the risk and occurrences of SA.

There were 168 positional candidate genes, representing 86 loci, associated with SA in cows bred by AI, of which 68 were previously identified as having a role in fertility (Table 3). Two positional candidate genes (FST5 on BTA17 and CDH13 on BTA18) have been identified as positional candidate genes associated with SA in a separate population of Holstein cattle and were associated with embryonic loss (3, 20). The cadherin 13 (CDH13) gene, identified as a positional candidate gene for SA in AI bred cows, encodes a calcium-dependent cell adhesion, cadherin protein that protects vascular endothelial cells from apoptosis due to oxidative stress (46, 47). CDH13 plays a role in cell migration and acts as a transcription repressor as a tumor suppressor gene (48–50). CDH13 is correlated with trophoblastic invasion anomalies that lead to fetal growth restriction in humans, which can cause fetal demise (51, 52). Other studies have found that CDH13 plays a role in early reproductive development and cholesterol biosynthesis (46). There were 23 positional candidate genes for SA in AI bred cattle that have previously been identified as enriched for embryonic loss and one (EYA2) that was enriched for SA in Holstein cattle (3, 20, 31). An additional 14 genes (SOX5, LOC101904842, LOC100847115, EYA2, LOC101905203, ZMYND8, LOC104973934, LOC112449409, CDH13, PHACTR1, AP1M1, LOC112447625, FAM32A, and TRNAE-UUC) were reported as positional candidate genes associated with SA in Holstein heifers bred by AI (13).

The four positional candidate genes of the most significant loci include the splicing factor proline and glutamine-rich (SFPQ) gene, which encodes an RNA-binding protein that is critical for maintaining transcriptional elongation of long genes and enabling histone deacetylase binding activity. Histone deacetylases are enzymes that catalyze the removal of acetyl functional groups from lysine residues of both nonhistone and histone proteins (53). The removal of epigenetic modifications regulates chromatin structure and transcription, while the deacetylation of non-histones controls a variety of cellular processes (53).

The second positional candidate gene from the most significant loci for SA is the RNA binding motif single-stranded interacting protein 3 (RBMS3) gene, which is involved in DNA replication, gene transcription, cell cycle progression, and apoptosis (54). This gene influences the regulation of microRNA expression or stabilization, inhibits the Wnt/β-catenin pathway, and affects other epithelial–mesenchymal transition-related transcription factors (55).

The third positional candidate gene is DLG associated protein 1 (DLGAP1), which is a part of the DLGAP family of proteins that act as scaffold proteins in the brain (56). These proteins play a vital role in synaptic scaling by regulating the turnover of both metabotropic and ionotropic glutamate receptors in response to synaptic activity (56). DLGAP1 is specifically involved in megakaryocyte biology and platelet function, potentially giving a proliferative advantage in hematopoietic cells (57).

Lastly, the fourth positional candidate gene is LOC107131937, a pseudogene of the small ribosomal subunit protein uS12 that is conserved across all domains of life (58). The ribosomal protein S12 is a component of the decoding center of the 30S ribosomal subunit and is involved in transfer RNA selection (59). These transfer RNAs are key in synthesizing proteins, serving as an adapter or link between the messenger RNA (mRNA) molecule and the growing chain of amino acids that make up the protein (60). Pseudogenes can be either functional or non-functional; the status of this pseudogene is currently unknown. If it is functional, it could impact the synthesis of important proteins that affect SA. Holstein infertility haplotypes have also been discovered, which are significant due to their association with infertility, embryonic loss, and stillbirth in Holstein cattle. This makes it important to determine if any SA-associated loci found in this study may also correspond to several identified Holstein infertility (HH) haplotypes. Several Holstein infertility haplotypes (HH) have been identified and are currently being tested to reduce infertility. These haplotypes include HH1 to HH6 (61) and are located on BTA5 at 63 Mb (HH1), on BTA1 at 94 Mb (HH2), on BTA8 at 95 Mb (HH3), on BTA1 at 1.3 Mb (HH4), on BTA9 at 92 Mb (HH5), and on BTA16 at 29 Mb (HH6) (62, 63). The nearest locus associated with SA in primiparous cows bred by AI was approximately 1.2 Mb away from HH5 on BTA 9. The next closest loci were four loci on BTA 8 that ranged from 8 Mb to 22 Mb away from HH3. Additionally, there was one locus that was 19 Mb away from HH2 on BTA 1, one locus that was 13 Mb away from HH6 on BTA16, and another that was 22 Mb away from HH1 on BTA5. Any other loci associated with SA were over 22 Mb away from these identified Holstein haplotypes, indicating independent regions.

Sixteen positional candidate genes were associated with SA in cows that were ET recipients, with five (SLC23A2, IDI1, PRNP, RASSF2, and GTPB4) previously identified as having roles in fertility. Solute carrier family 23 member 2 (SLC23A2) encodes a sodium-dependent vitamin C transporter linked to spontaneous preterm birth in women, with research suggesting that genetic variants may increase the risk of preterm delivery (64). IDI1 was previously reported as a leading-edge gene enriched with embryonic loss in Holstein heifers (20). Ras association domain family member 2 (RASSF2) is a tumor suppressor gene that regulates Ras signaling (65). The GTP binding protein 4 (GTPBP4) encodes a small GTP-binding protein located in the nucleus and is involved in signal transduction (66). There were no loci identified as associated with SA in primiparous cows receiving ET that were on the same chromosome as any of the six Holstein infertility haplotypes mentioned earlier.

Three genes (RASSF2, GTPBP4, and IDI1) have been previously reported to be differentially expressed in the endometrium of pregnant beef heifers compared to non-pregnant, high-fertility beef heifers (32). IDI1, PRNP, and LARP4B were expressed in the reproductive tissues of cattle with pregnancies generated via embryo transfer, thus supporting their role in the function of a normal pregnancy (33).

Research is ongoing to determine the genomic differences between SA in AI-bred and ET-recipient cattle. There was no overlap of loci or genes associated with SA in cows bred by AI with those associated with SA in ET-recipient cows. This lack of overlap could be due to insufficient power (power = 0.2) to detect an association or due to distinct physiological, biological, and genomic regions essential for pregnancy establishment and maintenance in AI-bred and ET recipients. Lower statistical power also increases the risk of false negatives in the analysis; therefore, further studies need to be conducted with larger sample sizes to confirm these results. To the author's knowledge, there are no other association studies that have investigated loci associated with fetal loss in ET-recipient cattle, but there have been studies investigating embryonic loss in ET recipients. Similar results have been identified in loci associated with embryonic loss in cows and heifers bred by AI and those that were ET recipients (67) and are consistent with physiological studies that have identified differences between AI-bred and ET-recipient pregnancies.

For primiparous cows bred by AI and those receiving ET, the greatest number of SA (as indicated by cows returning to estrus) occurred on day 45 of gestation. This could indicate possible issues or deformities of the placenta that lead to SA. Functional insufficiency of the placenta often results in increased offspring mortality (68). The gestational interval between days 28 and 60 is a pivotal period for the development of the placenta (69). This period of gestation requires substantial growth and differentiation of embryonic and trophoblastic cells, with significant communication between the dam and the embryo occurring (69). After day 30 of gestation, a firm connection between fetal cotyledons and maternal caruncles occurs, with the fetus developing a distinguishable adult organ system by day 60. Fetal loss during this period can result from developmental deficiencies such as inadequate placentation, difficulty transitioning from amniotic to allantoic nutrition, changes in placental vascularization, and underdevelopment of the embryo or fetus (69). Twenty-five percent of pregnancies from in vitro-produced embryos are reported to be lost between days 28 and 60 of gestation (69). In these embryos, there is a clear underlying physiology of placental abnormalities that contribute to the high incidence of pregnancy loss during this period.

Other factors can impact SA, such as infectious diseases, nutritional status, and overall management. In this study, we accounted for infectious diseases and some nutritional issues such as ketosis by excluding animals from the study that were recorded as experiencing these events during their pregnancy. However, not everything can be carefully managed or recorded in the dairy's records (such as body condition score), so there is potential for confounding variables to go undetected and potentially impact the occurrence of SA.

A future direction for the analysis of SA with a greater number of cases would be to split the gestational window into two categories to include early losses (approximately 42 to 120 days) and late losses (121–260 days). As reported in the results, there was a greater number of SA that occurred between days 42 and 120 for both primiparous cows bred by AI and those that were ET recipients. Further analysis focusing on SA occurring between days 42 and 120 could also be beneficial, as early losses are more often linked to genetic factors than late-term losses. A future study examining these early and late losses could be very helpful in identifying genetic regions and biological components that impact SA losses at different times during gestation.