Statistical methods to disentangle genetic effects influencing infertility and early fetal viability with a genome-wide application.

Recent estimates from the World Health Organization suggest that 17.5% of the world’s population experience infertility at some point during their lifetime [ 1 , 2 ]. Infertility is usually defined as the inability to achieve a clinically recognized pregnancy after having tried for 12 months or more [ 3 ]. Even under optimal conditions, the maximum chance of achieving a clinically recognized pregnancy in a given cycle is around 30-40% [ 4 ]. Given that an increasing number of couples now choose to delay childbearing until later in life [ 5 , 6 ], which may negatively impact their reproductive success, it is not surprising that the use of assisted reproductive technology (ART) has also increased in many parts of the world [ 7 ]. Recent data from the Medical Birth Registry of Norway show that approximately 5% of all newborns in 2022 were conceived by ART [ 8 ]. Genetic effects influencing infertility may operate through a variety of biological mechanisms, including chromosomal abnormalities, monogenic disorders, or traits with complex inheritance [ 9 – 11 ]. Several studies have reported associations between single-nucleotide polymorphisms (SNPs) and infertility-related conditions, such as endometriosis and ovulation disorders in women and aberrant sperm characteristics in men [ 12 – 16 ], while others have screened for genetic variants associated with unexplained (idiopathic) infertility [ 17 , 18 ]. A few recent studies have highlighted the importance of considering recessive models to search for genetic effects linked with reproductive traits. For example, Ruotsalainen et al. identified SNPs in the TBPL2 gene with genome-wide significant recessive effects on female infertility [ 19 ]. Oddsson et al. discovered that women homozygous for a stop-gain (nonsense) variant in the CCDC201 gene reached menopause on average nine years earlier compared to other women [ 20 ]. Most of the above studies have examined male and female reproductive traits separately. While they, along with other research, have provided significant insights into the genetic determinants of infertility, the cause remains unknown in approximately 15-30% of couples struggling to conceive [ 21 – 23 ]. In some of these couples, it is plausible that specific combinations of male and female germ cells lower the likelihood of fertilization. A couple’s use of ART to conceive may be regarded as a proxy for their infertility. In an earlier study, we compared ART and non-ART pregnancies to identify male-female allelic interaction effects influencing a couple’s ability to conceive beyond what would be expected from each partner’s individual risk of infertility, revealing interaction effects with SNPs in DNAH17 , a gene with known links to infertility [ 24 ]. Even after a successful fertilization, a substantial proportion of embryos may still perish at an early stage before the pregnancy is clinically detectable. Around 60% of all conceptuses have been estimated to be lost prior to the missed menstrual period [ 4 , 25 ]. Repeated early pregnancy losses may also be perceived as couple-level infertility, and parental alleles associated with infertility may in reality increase the risk of early pregnancy loss after being transmitted to the embryo. Currently, research on genetic determinants of embryonic survival is limited. However, in a study by Mayor-Olea et al., it was suggested that selection in favor of the T allele in the MTHFR 677C>T polymorphism may be associated with early fetal viability [ 26 ]. Recently, Arnadottir et al. provided more direct evidence that at least 1 in 136 pregnancies was lost due to a pathogenic small sequence variant, based on whole-genome sequencing of samples from fetal-parent triads affected by pregnancy loss and adult control triads [ 27 ]. To shed further light on infertility caused by early pregnancy losses, we previously investigated allele distributions in ART-conceived offspring, conditional on parental genotypes [ 28 ], and analyzed case-parent units using the R package Haplin [ 29 – 31 ]. If genetic influences on early pregnancy loss are more pronounced among couples who seek ART treatment, it is conceivable that the ART procedures mitigate such genetic influences. Thus, genetic variants associated with early pregnancy loss would exhibit allele frequencies that deviate from expected Mendelian parent-to-offspring transmission patterns in both non-ART and ART offspring. Under such a scenario, the frequencies of potentially deleterious alleles would also be expected to differ between non-ART and ART offspring. Indeed, we found several alleles whose frequencies in the ART offspring differed from expected proportions in two genes, CXXC4-AS1 and DYNLRB2-AS1 [ 28 ]. A research question that remains to be explored is how to identify whether specific genetic variants influence the probability of fertilization or embryonic survival. When parental and fetal genetic effects are analyzed separately, each type of effect may be “aliased” by the other; in other words, a parental effect could manifest as a fetal effect, and vice versa [ 32 ]. To address this analytic problem, we developed a novel approach to disentangle parental genetic interaction effects influencing fertilization success from fetal genetic effects influencing early fetal survival. Our methodology is based on log-linear modeling and builds on existing framework for case-parent triad analyses available in Haplin [ 29 ]. We conducted genome-wide analyses of the fetal effect on fetal survival, with and without adjustment for parental interaction effects. An important note on terminology is that our dataset contains genotype data only for offspring who reached a gestational age well beyond the embryonic stage. Therefore, while our theoretical framework is centered around “embryo survival”, our analyses focus on “fetal survival”, and we use both terms accordingly. More details on data inclusion are provided later. Similarly to how fetal effects may differ between non-ART and ART offspring, parental effects may also differ in non-ART compared to ART parents. Therefore, separate analyses were performed in two study samples of case-parent units: one with naturally conceived offspring and the other with offspring conceived through ART. We provide details on the model implementations, along with the results of the genome-wide analyses.

Fig 3 shows the results of the analysis of single-dose fetal effects on fetal survival in the non-ART sample, without adjustment for parental interaction effects. Results for the corresponding double-dose fetal effect are presented in S1 Fig . The Manhattan plot in Fig 3 A shows one distinct peak of neighboring SNPs on chromosome 6. The most significant SNP, rs114918746, had a p -value of 4.44 × 10 −16 , three additional SNPs (rs113299093, rs2517588 and rs2905736) had p -values below the Bonferroni-corrected genome-wide significance threshold of 5 × 10 −8 , and a fifth SNP (rs114011376) had a p -value just slightly above the genome-wide threshold. Applying an FDR threshold of 0.1, we identified 10 additional significant SNPs, bringing the total to 15, including the five lead SNPs. The RR and MAF estimates for these SNPs are provided in Table 1 . Most of the estimated MAFs for these SNPs were around 1-3%. However, a few more common variants were identified on chromosomes 12, 16, and 17, and their MAFs ranged from 5.6% to 11.5%. All RR estimates were below 1, ranging from 0.75 to 0.92. As indicated by the quantile-quantile (QQ) plot in Fig 3 B, the p -value distribution was negligibly inflated ( λ = 1.107). Also included in Table 1 are the HWE p -values, and there was considerable variation in these p -values across the 15 significant SNPs, ranging from 0.006 to 0.993. The smaller p -values indicate violations of the HWE assumption. At a threshold of 0.05, three of the five lead SNPs (rs2517588, rs113299093, and rs114918746) deviated from the expected HWE distribution. The lead SNPs are located in or near genes such as Fibroblast Growth Factor Receptor 3 Pseudogene 1 ( FGFR3P1 ), Ribosomal Protein L7 Pseudogene 4 ( RPL7P4 ), Mediator of DNA Damage Checkpoint 1 ( MDC1 ), HLA Complex Group 17 ( HCG17 ), HLA Complex P5 ( HCP5 ), MHC Class I Polypeptide-Related Sequence B ( MICB ), and Notch Receptor 4 ( NOTCH4 ). Panel A displays the negative log 10 -transformed p -values and the red line indicates the Bonferroni-corrected genome-wide significance threshold of 5 × 10 −8 . Panel B shows the corresponding quantile-quantile (QQ) plot, along with the genomic inflation factor ( λ ). 1 Only one SNP from each pair in strong linkage disequilibrium ( r 2 >0.95) is presented in this table. 2 The q -values are estimated from the observed p -values using the Benjamini-Hochberg FDR method. 3 Hardy-Weinberg equilibrium test p -values. The same SNPs presented in Table 1 were re-analyzed in the non-ART sample, but this time adjusting for parental complementary and parental threshold-based interaction effects. The fetal effect RR estimates changed negligibly after these adjustments (top panel of Fig 4 ). The same SNPs were subsequently analyzed in the ART sample, first without adjustment for parental interaction effects and then with adjustment for parental complementary and parental threshold-based interaction effects (bottom panel of Fig 4 ). Most of the identified fetal effects were no longer statistically significant when analyzed in the ART sample. However, one SNP, rs113299093, showed a significant unadjusted fetal effect in the ART sample, similar to that observed in the non-ART sample. Its p -value of 1.09 × 10 −3 was below the candidate-SNP significance threshold of 0.05/15, based on the re-analysis of 15 SNPs. The RR estimates and p -values for the fetal and parental interaction effects of the SNPs listed in Table 1 , analyzed in the models described above for both the non-ART and ART samples, are provided in S1 Table . A fetal effect RR estimate greater than 1 indicates an increased “risk” of fetal survival for carriers of the minor allele, whereas interpretation of parental interaction RRs depends on whether they are estimated within the ART or non-ART sample. In the ART sample, an RR above 1 corresponds to an increased risk of needing ART to conceive, i.e., a potentially reduced reproductive capacity, for couples carrying the interaction genotype combinations. Conversely, in the non-ART sample, an RR above 1 indicates an increased risk of not needing ART, reflecting enhanced reproductive capacity. Manhattan and QQ plots for the genome-wide analyses of parental interaction effects in the non-ART and ART samples are presented in Figs 5 and 6 . Corresponding plots for the single and double-dose fetal effects from the same models are presented in S2 – S6 Figs . While not genome-wide significant, some noteworthy clusters of neighboring SNPs were apparent in these figures. Interestingly, these results also partially overlapped with those in our previous work [ 28 ]. More specifically, S5 A and S6A Figs are comparable to Fig 2A in [ 28 ]. For the six SNPs in Table 1 that were not included here, the parental interaction effect RRs in the ART sample could not be estimated properly, likely due to model convergence issues as a result of the small sample size. Panels A and C display the negative log 10 -transformed p -values, and the red line indicates the Bonferroni-corrected genome-wide significance threshold of 5 × 10 − 8 . Panels B and D show the corresponding quantile-quantile (QQ) plots, along with the genomic inflation factor ( λ ). Panel A shows the negative log 10 -transformed p -values, with the red lines indicating the Bonferroni-corrected significance threshold of 5 × 10 −8 . Panels B and C show the corresponding quantile-quantile plots, along with the genomic inflation factor ( λ ). Although some caution is warranted when interpreting the validity and significance of the parental interaction effect estimates, we nevertheless opted to present these results here to motivate a methodological discussion on modeling strategies for such genetic effects. The analyses of parental complementary and parental threshold-based interaction effects in the non-ART sample identified 32 and 27 genome-wide significant SNPs, respectively; however, none of these SNPs overlapped. Details of these results, including MAF and RR estimates, as well as p -values, are provided in S2 Table . Most of the identified significant SNPs were rare, and had RR estimates mostly ranging from 1 to 3. In this context, an RR greater than 1 indicates that the specific parental genotype combination under study is more common than expected among parents who conceived without the use of ART. The significant SNPs associated with complementary interaction effects were broadly scattered across the genome, whereas some of the SNPs significantly associated with threshold-based interaction effects were more tightly clustered, particularly on chromosome 16. These latter SNPs are located in two genes, Zinc Finger CCHC-Type Containing 14 ( ZCCHC14 ) and Junctophilin 3 ( JPH3 ), neither of which have any previously published connections to infertility. Among the remaining significant parental interaction effects in the non-ART sample, two SNPs (rs35974259 and rs360673) are particularly noteworthy in that they are located in or near genes relevant to infertility. SNP rs35974259 on chromosome 13 is associated with a parental complementary interaction effect and is located near the gene for Sperm Acrosome Associated 7 ( SPACA7 ). The protein product of this gene assists spermatozoa in breaking through the protective coating of cells surrounding the targeted oocyte during fertilization [ 46 ]. SNP rs360673 on chromosome 1 is associated with a parental threshold-based interaction effect and resides within Synaptonemal Complex Protein 1 ( SYCP1 ), a gene known to play a crucial role in meiotic chromosome pairing and recombination [ 47 ]. As shown in the Manhattan and Miami plots in Figs 5 and 6 , multiple SNPs were significantly associated with parental complementary and/or parental threshold-based interaction effects in both the non-ART and ART samples. There was no overlap between the two samples, as shown in the Miami plots directly comparing them in S8 Fig . The accompanying QQ plots revealed potential issues with the models; notably, the p -value distributions from the analyses in the non-ART sample were inflated, with genomic inflation factors of λ = 1.383 and λ = 1.364, respectively. On the other hand, for the results in the ART sample, while visual inspection of the QQ plots also indicated varying degrees of inflation, the corresponding λ -values were both less than 1 (0.697 and 0.828, respectively), suggesting that the overall distribution was deflated . This was due to a disproportionately large number of p -values near 1, which was not evident from a visual inspection of the QQ plot. To explore these issues further, SNPs with a MAF below 5% were excluded and new QQ plots were generated (see S9 Fig ). Excluding the rarest variants resulted in the inflation factors for the ART analyses being closer to 1 (1.093 and 1.099). Nevertheless, the results of the non-ART analyses remained substantially inflated, suggesting that the problem with the model may be due to factors other than low MAFs. Within the ART sample, the results of the parental complementary and parental threshold-based interaction effects showed a greater degree of overlap, as illustrated by the Miami plot in Fig 6 . Notably, one cluster of top SNPs on chromosome 7 was shared by both types of parental interaction effects. Further details, including RR estimates and p -values, for the SNPs with genome-wide significant parental interaction effects in the ART sample—either complementary, threshold-based, or both—are presented in S3 Table . Most of the RR estimates were quite large, which is likely due to low MAFs and small numbers of parental couples carrying the specific allele combinations. Nevertheless, RR estimates greater than 1 again indicate increased frequencies of the specific genotype combinations relative to those expected within this sample, i.e., they are more common among parents who used ART to conceive. Further investigation of these results revealed several relevant genes in the regions encompassing these SNPs, as shown in the regional association plots in Fig 7 and S7 , including Actin Beta ( ACTB ), Fascin Actin-Bundling Protein 1 ( FSCN1 ), and Ring Finger Protein 216 ( RNF216 ). The green arrows indicate the direction of transcription for the genes in the region, while recombination rates are represented by blue lines. The color gradient reflects the degree of linkage disequilibrium, with associated r 2 values. We also conducted simulation-based power calculations for the model with adjustment for parental complementary interaction effects. Based on the framework of the hapSim function [ 48 ] in the Haplin R package, we simulated offspring-parent SNP data under scenarios with true parental complementary interaction effects ranging from RR = 1.2 to RR = 2.2, and with the fetal effect fixed at the null (RR = 1). We then analyzed the simulated data and estimated statistical power as the proportion of analyses in which a parental complementary interaction effect was significant at the Bonferroni-corrected threshold of 5 × 10 −8 . Varying both MAF and sample size, we found that in a sample with a size comparable to our ART sample, the estimated power to detect a SNP with MAF = 10% and RR = 2.2 was 22%. If the MAF increased to 0.15, the estimated power increased to 73%. In a larger sample, with a size comparable to our non-ART sample, the estimated power to detect a SNP with MAF = 5% and RR = 1.5 was 79%. Results from these power calculations are presented in S10 Fig .

In this study, we investigated parental and fetal genetic effects influencing infertility and fetal survival in family triads and dyads from the MoBa cohort. Case-parent units were categorized according to whether the offspring was conceived using ART or not, and separate analyses were conducted in these two groups. In the non-ART sample, we identified genome-wide significant fetal effects on fetal survival, primarily for several variants located on chromosome 6. These associations remained significant after adjustments for parental complementary and parental threshold-based interaction effects, confirming that these effects are fetal in origin. When re-analyzed in the ART sample, the results were either consistent with those in the non-ART sample or inconclusive. In other words, there was not enough evidence to conclude that the fetal effects observed in the non-ART offspring are absent in the ART offspring. One possible interpretation is that the observed fetal effect on fetal survival is not directly influenced by ART use. Many of the identified SNPs had low MAFs, which is as expected, as alleles with potentially detrimental effects on reproductive ability are unlikely to be common in the population. The RR estimates were all below 1, indicating that the effect (minor) alleles in the offspring were even less frequent than expected given the parental genotypes and Mendelian inheritance. The potential negative selection of these alleles suggests that they are associated with a lower probability of fetal survival. Furthermore, three of the SNPs with a genome-wide significant fetal effect deviated from HWE. A possible explanation for this is that the distribution of alleles associated with reproductive ability may shift over generations. Overall, due to the family-based study design, the fetal effect is inherently robust against population stratification, supporting the validity of these results and minimizing the risk of bias, for instance, from assortative mating. Among the genetic variants identified as being associated with a fetal effect on fetal survival, several are located in or near genes involved in key reproductive processes, including fertility regulation, immune mechanisms, ovarian function, and early pregnancy development. For example, MDC1 is linked to male fertility through its regulation of meiotic progression in spermatocytes [ 49 , 50 ]. MICB has been suggested to play a role in immune-related fertility mechanisms [ 51 ]. HCP5 , a key long non-coding RNA (lncRNA), is implicated in premature ovarian insufficiency and regulates critical pathways involved in ovarian function and fertility [ 52 ]. NOTCH4 is essential for several reproductive processes, including embryo implantation, placentation, and uterine function [ 53 , 54 ]. These genes are therefore relevant not only for the outcome under study—fetal survival—but also for other related phenotypes, such as the parents’ ability to conceive. One possible explanation for the observed fetal effects in genes relevant for parental phenotypes is that alleles increasing the risk of infertility in parents may also, if passed down to the offspring, influence the probability of fetal survival. Furthermore, although the fetal effects were not attenuated after adjusting for parental complementary and parental threshold-based interaction effects, we cannot exclude the possibility that other types of parental effects could be involved. We queried the 15 SNPs from Table 1 in the online database dbSNP [ 55 ] and found that several are classified as intronic variants, while others have no assigned functional consequence. Intronic variants may still affect gene expression by influencing mRNA stability, polyadenylation, miRNA binding, or other post-transcriptional regulatory processes. However, detailed functional and epigenetic annotations remain limited for many of these variants, and further investigation is needed to clarify their potential biological impact. The parental interaction effects, as implemented in the models described here, are assumed to influence the overall ability to achieve fertilization, i.e., a study outcome that differs from fetal survival. When analyzed within the same model, the aim is to separate fetal and parental effects. While the fetal effects capture the direct influence of the fetal genotype on fetal survival, the parental interaction effects, adjusted for fetal effects, reflect the couple’s ability to achieve fertilization. For example, if a parental complementary interaction effect is identified with an RR estimate greater than 1 in the ART sample, the interpretation is that parents carrying the specified complementary genotype combination occur more frequently among ART triads than expected. In other words, such a result suggests an association between the parental complementary genotype combination and reduced reproductive capacity, assuming random mating. Since the parental interaction effect estimated in the non-ART sample reflects the opposite outcome—namely, enhanced reproductive capacity—one might expect the parental interaction RR estimates for the same SNPs to show opposite directions in the ART and non-ART samples. As shown in Fig 4 and S1 Table , this pattern occurred occasionally but was not consistent. Furthermore, it can be hypothesized that genetic variants detrimental to parental ability to achieve fertilization may also negatively impact the likelihood of fetal survival. If so, fetal and parental interaction effect estimates would be expected to align in direction within the non-ART sample, but diverge in the ART sample. However, this expectation was not consistently supported by the results presented in Fig 4 and S1 Table . Compared with the fetal effect, parental interaction effect estimates are more susceptible to bias in case-triad analyses. Essentially, the four assumptions stated for the model definition in Sect 2.3 are not all necessarily met. The first assumption, that mothers and fathers originate from populations with equal allele frequencies, is not directly tested with our data. However, principal component analysis was performed on the full MoBa genetic dataset by Corfield et al., establishing that approximately 95% of MoBa participants are of European descent [ 41 ]. Nevertheless, MoBa participants represent a wide geographic range within Norway, and distinct geographic genetic subgroups have been identified in the Norwegian population [ 56 ]. In summary, we cannot claim that the first assumption is fully met, but any violations are likely limited in extent or effect. The second model assumption concerns HWE, which we tested directly in our data. Because the HWE test pertains to the underlying data structure, it is only minimally influenced by which parameters are estimated. In other words, the HWE test p-values presented for the non-ART fetal effect results in Table 1 can also be used to assess the validity of the non-ART parental interaction effect estimates for the same SNPs. A few of the SNPs in Table 1 exhibited signs of departure from HWE, which could introduce bias into the corresponding parental interaction estimates. Regarding the third model assumption of Mendelian transmission, it is important to recognize that deviations from the expected random segregation of alleles can occur at several stages, including fertilization or earlier. In our current analyses, such deviations may manifest as apparent genetic associations with early fetal viability. Wilcox et al. highlighted the assumption of Mendelian inheritance in case-parent triad studies, emphasizing its importance in accurately identifying genetic variants associated with disease in the offspring [ 57 ]. In the current study, we demonstrated that this assumption may not always hold and should therefore be specifically addressed in case-triad analyses. The fourth model assumption listed in Sect 2.3 , relates to random mating. This condition is commonly hard to fulfill in genetic association studies such as ours, and previous research has indeed demonstrated assortative mating in several traits for couples in the MoBa cohort in particular [ 58 ]. While such violations are less problematic for the estimation of fetal effects, they may introduce bias in the parental interaction effect estimates. Moreover, the genetic selection of parents into the non-ART or ART groups may depend on factors unrelated to infertility, potentially biasing the results. If specific genotype combinations accumulate in couples as a result of positive or negative assortative mating, population stratification, or selection bias, they may manifest as false-positive parental interaction effects. This is likely a greater issue for SNPs with low MAFs, as the interaction depends on two loci rather than one, and genotype frequencies for these SNPs have the potential to deviate more relative to the expected frequencies under Mendelian inheritance. Nevertheless, variants that influence parental ability to conceive, i.e., their fecundity, are likely to be rare in the population. The analyses of parental complementary and parental threshold-based interaction effects primarily identified lead SNPs with low MAFs in both the large non-ART sample and the substantially smaller ART sample. The QQ plots for the parental interaction effects in the non-ART sample revealed model-related issues that could not be resolved by excluding low-frequency variants. Interestingly, the impact of excluding rare variants differed between the non-ART and ART samples, with QQ plots showing improved model calibration after excluding SNPs with low MAF in the ART sample only. Since this effect was not observed in the non-ART sample, it could simply have been due to chance. A possible explanation could be that the impact on genomic inflation caused by confounding factors could be larger in a larger sample size. Analyses of parental interaction effects in the ART sample are likely subject to similar overall issues as those in the non-ART sample. Furthermore, certain subtypes of infertility may be polygenic, which could contribute to the genomic inflation observed in the QQ plots. Polygenicity could be assessed using methods such as LD score regression [ 59 ]. However, when the unit of analysis is parental couples, as in this study, no immediately applicable LD score regression approach is available. Although the interpretation of the parental interaction analyses is not straightforward, the choice of parameterizations still merits further discussion. The complementary and threshold-based parameterizations presented here are just two examples of how parental incompatibility may manifest as a genetic effect affecting reproductive ability. Several other parameterizations could also be used to model relevant parental interaction effects, some of which were explored in case-control analyses in our previous work [ 24 ]. There, we showed that different parameterizations can yield distinct results, underscoring the value of a broad analytical approach incorporating multiple models. The choice of model may be guided by prior hypotheses about the underlying biological mechanisms. For instance, if a parental genetic incompatibility is suspected to influence reproductive capacity, the complementary model may be particularly appropriate. Overall, the complementary and threshold-based parameterizations of the parental interaction effect exhibit some overlap. Parents who harbor the complementary interaction genotype combinations—where one parent is homozygous for the effect allele and the other for the reference allele—also meet the criteria for the threshold-based interaction combinations, where the combined number of effect alleles in both parents is required to be 2 or higher. In other words, a true complementary effect could be detected using a threshold-based parameterization, and vice versa. This may not be true for other threshold values. For example, a threshold of 3 excludes the complementary genotype combinations, while a threshold of 1 includes all parental genotype combinations except the combination where both parents are homozygous for the reference allele. A promising direction for future research could be to investigate when and why the complementary and threshold-based models yield overlapping results—and when they do not—as this could help inform model choice in other settings. It is also important to acknowledge that additional genetic mechanisms not captured by our models may influence the outcomes. In the context of infertility and fetal survival, a maternal-fetal genotype incompatibility, as described by Sinsheimer et al. may also be relevant [ 60 ]. Models accounting for maternal-fetal interactions, as well as other types of genetic effects, have been summarized by Ainsworth et al. [ 61 ]. Typically, attempts to estimate multiple genetic parameters are constrained by the available degrees of freedom. Overlapping results between the two parental interaction models were indeed observed in the ART sample, where several of the same significant SNPs were identified as having both parental complementary and parental threshold-based interaction effects. Of note, these SNPs are located in a region containing several genes associated with infertility. For instance, ACTB is involved in proper sperm formation and male fertility through its role in organizing the mitochondrial sheath during spermiogenesis [ 62 , 63 ]. FSCN1 , an actin filament-bundling protein, is also essential for male fertility and spermatogenesis [ 64 , 65 ]. Similarly, RNF216 encodes an E3 ubiquitin ligase that plays a crucial role in male fertility through regulating spermatogenesis, meiosis, and PKA stability in the testes [ 66 , 67 ]. The phenotypic relevance of these genes supports the notion that these associations are true rather than spurious. For the non-ART analyses, no significant SNPs overlapped between the complementary and threshold-based models. One possible explanation for why an overlap in results for the parental interaction effect was observed only in the ART sample is that the parental phenotype studied in this group—risk of needing to use ART—is more specific and well-defined than the contrasting phenotype examined in the non-ART sample. Another possibility is that if some of the findings are influenced by confounding factors, such as assortative mating, these factors may be more pronounced in the ART group. The high quality of the genetic and registry data combined with the large sample size, particularly for the non-ART sample, are among key strengths of this study. There are likely only a few other datasets in the world where a genome-wide case-parent triad analysis at this scale is feasible. In contrast, although the ART sample also ranks among the largest such datasets globally, the sample size remains relatively modest in a genome-wide context, limiting the statistical power to detect small or moderate effects in this group. The substantial difference in sample size between the ART and non-ART samples limits the reliability of comparisons between these groups. Furthermore, the relatively low initial participation rate in MoBa suggests that the cohort may not be fully representative of the general population. For instance, MoBa participants generally have higher levels of education compared to the broader population [ 68 ]. However, Nilsen et al. analyzed MoBa data and demonstrated that bias in prevalence estimates does not necessarily imply bias in the estimated associations between exposures and outcomes [ 69 ]. During the years of recruitment into the MoBa cohort (1999–2008), the overall Norwegian population grew [ 8 ], and ART use increased modestly both in the population [ 8 ] and within the cohort (See S11 Fig ). However, the impact on our analyses is likely minimal given the relatively homogeneous ancestry of the cohort. As is often the case, genotypes from pregnancy losses were unavailable and therefore could not serve as a suitable control group for these analyses. Given this limitation, the case-parent triad design utilized here represents one of the most robust approaches available. Potential analyses to follow up on or complement our work include studies using other measures of reproductive ability, such as time to pregnancy, and functional analyses of the genetic associations reported in this study. In conclusion, we identified fetal effects on fetal survival that persisted after adjusting for parental interaction effects in the non-ART sample. The evidence was inconclusive regarding whether these effects were also present in the ART sample. Nevertheless, the reliability of the observed fetal effects in the non-ART offspring is supported by the sample size, the robustness of the fetal effect estimates due to the case-triad study design, and the relevant genetic context. We also observed parental interaction effects in previously reported infertility-related genes, especially in the ART sample. As with most observational studies, our findings may have been influenced by biases that are difficult to account for in the models. This study is among the first to examine the joint estimation of parental interaction effects and fetal effects using a family-based study design. The identification of numerous genes related to reproductive ability supports the validity of the approach. Future studies are nevertheless needed for replication and validation of our findings.

The study was conducted using data from the Norwegian Mother, Father, and Child Cohort Study (MoBa). The establishment of MoBa and initial data collection was based on a license from the Norwegian Data Protection Agency and approval from The Regional Committees for Medical and Health Research Ethics. The MoBa cohort is currently regulated by the Norwegian Health Registry Act. The current study was approved by The Regional Committees for Medical and Health Research Ethics (#2017/1362). The adult MoBa participants provided written statements of consent at inclusion, both for themselves and on behalf of their children [ 33 ] For a given couple, there may be multiple scenarios leading to a successful conception and having a baby, or trying without succeeding. To guide the statistical modeling, Fig 1 outlines a number of different reproduction trajectories that may lead to one of these two outcomes. The trajectory to reach a birth at a viable gestational age without the use of ART is marked with light gray boxes in Fig 1 and consists of the following three steps: 1) Each of the parents is independently biologically fertile; 2) as a couple, they are able to achieve fertilization without ART; and 3) the embryo survives, develops, and reaches maturity. Each of these steps can be disrupted by detrimental genetic or environmental factors. While environmental factors can, in principle, be adjusted for at each level, our modeling approach emphasizes how parental and fetal genes collectively contribute to the outcome. Each numbered stage is further illustrated with an example in Fig 2 . Step 1 is the part most commonly associated with infertility, as numerous genetic factors are known to lead to female or male infertility. Couples where either or both partners are affected by genetically influenced infertility may decide to use ART. Consequently, if the ART procedure is successful, and the pregnancy reaches viability, the couple will end up in our second group, consisting of ART births at a viable gestational age. The trajectory toward a viable birth after ART use is marked with light purple boxes in Fig 1 . Step 2, however, requires a certain level of genetic compatibility between the parents. Notably, while each parent might be able to conceive with other partners, the combined genetic makeup of both parents may increase the difficulty of fertilization beyond what would be expected from their individual, genetically determined reproductive potential alone. Previously, our group explored various models of compatibility, or interactions, between parental alleles [ 24 ], inspired by an article that summarized all possible two-locus, two-allele penetrance models for dichotomous traits [ 34 ]. As an example, the “complementary” model assumes that opposite parental genotypes, i.e., a maternal aa combined with a paternal AA , or vice versa, can impact the likelihood of fertilization. Step 3 primarily relates to embryo survival and early pregnancy losses. If a couple experiences frequent early pregnancy losses, they may decide to pursue ART treatment. An ART pregnancy could then be successfully established, progressing beyond the risk time window of early miscarriage, either because the ART procedure enhances early fetal survival or because the fertilized egg does not carry the deleterious allele combination responsible for previous early losses. While Steps 1 and 2 above relate to the parental genotypes, Step 3 mainly focuses on the fetal genotype. However, parental and fetal effects may be aliased, meaning that parental effects could be mistaken for fetal effects, and vice versa. Examples of these scenarios are illustrated in Fig 2 for each of the three steps. Solid cells represent the direct genetic effect, while patterned cells (diagonal lines) indicate genotypes whose frequencies may be altered by the direct genetic effects. For Step 1, consider a scenario where maternal infertility is associated with an elevated frequency of a recessive allele a at a specific SNP. In ART triads, maternal aa genotypes would then be expected to occur more frequently than paternal aa genotypes within the same triads, and also more frequently than maternal aa genotypes in non-ART triads. Moreover, the distribution of the fetal genotypes would also be shifted accordingly in both ART and non-ART triads. For Step 1 illustrated in Fig 2 , filled cells represent elevated maternal genotype frequencies in ART triads, resulting directly from genetic selection. Patterned cells indicate altered fetal genotype frequencies, which occur as an indirect consequence of this selection. For Step 2, Fig 2 illustrates the selection under the complementary model, where opposite combinations of paternal alleles (filled cells of the same color) are more prevalent among ART triads than expected, and are subsequently less prevalent in non-ART triads. A maternal genotype of AA paired with a paternal genotype of aa , or vice versa, will always result in a heterozygous child. Consequently, this complementary interaction effect is indirectly reflected in the altered frequencies of heterozygous genotypes among children (patterned cells). For Step 3, a deviation in the fetal allele distribution from what would be expected given the parental genotypes and Mendelian inheritance may suggest that the presence of a specific allele influences the probability of embryo survival. This phenomenon may occur in both ART and non-ART triads. However, couples experiencing frequent early losses due to the fetal genotypes being less compatible with early viability may proceed to ART treatment, and in such cases, the parents are more likely to carry the deleterious alleles themselves. ART triads may therefore exhibit distortions in both parental and fetal allele frequencies. For Step 3 in Fig 2 , this scenario is exemplified by a fetal recessive genetic effect, which directly impacts the fetal aa frequencies (filled cells), whereas the corresponding parental allele frequencies are indirectly impacted by the selection of couples into the ART group (patterned cells). From the above discussion, it is evident that ART triads reflect a mixture of two related study outcomes: 1) the risk of ART use at the parental level; and 2) survival at the embryo level. Previously, we modeled parental interactions in a parental case-control setting, defining parents who used ART as cases and those who conceived naturally as controls [ 24 ]. We also demonstrated how our R package Haplin [ 29 ] can be used to analyze fetal survival in ART triads, estimating fetal effects as well as maternal and parent-of-origin effects, but without modeling the simultaneous effect of parental interactions [ 28 ]. In the current study, we extended the Haplin methodology to account for possible overlaps between parental and fetal genetic effects by jointly modeling parental interactions alongside fetal effects. By focusing the analyses on ART triads and non-ART triads separately, the risk of bias from population stratification is reduced, which allows for a more precise model. We will show that while within-triad analyses may still be exposed to biases from population substructure, such biases are less likely to have a serious impact on the results. At a specific multi-allelic locus, assume that there are H different alleles present in the population, and let A 1 , A 2 , … , A H denote those alleles. Let A i A j be the genotype of the mother in the triad, and A k A l the genotype of the father. Here, i , j , k , l ∈ { 1 , 2 , … , H } , and for instance i  =  j  = 2 would mean that the mother is homozygous with the A 2 A 2 genotype. In the mother-father-child triad context we need to distinguish between those alleles that are transmitted to the child and those that are not. We write M = ( A i , A j ) and F = ( A k , A l ) for the ordered genotypes of the mother and father, respectively, where we assume that the first allele from each parent is not transmitted, whereas the second allele is transmitted to the offspring; hence, the offspring genotype can be expressed as C = ( A j , A l ) . To describe the population genotype distributions, we will for instance write P ( M ) = P ( ( A i , A j ) ) for the proportion of mothers having the (ordered) M genotype, i.e., as a simplified notation for the probability that a given mother has the ordered genotype M , and similarly P ( M , F ) for the proportion of triads where the mother has the M genotype and the father the F genotype. The distribution of parental alleles in the background population, P ( M , F ), is particularly important since it serves as a baseline for comparison when looking at the distribution actually observed in, for instance, the ART triads. For the modeling, we will consider the consequences of the following basic assumptions: 1) Mothers and fathers originate from populations with equal allele frequencies, with p 1 , p 2 , … , p H denoting the (parental) background population frequencies of the A 1 , A 2 , … , A H alleles, and thus ∑ h = 1 H p h = 1 ; 2) The parental population is in Hardy-Weinberg equilibrium (HWE), which implies that the unordered maternal genotype A i A j has population frequency 2 p i p j when i ≠ j and p i 2 when i  =  j , and similarly for the paternal genotype; 3) There is Mendelian transmission of alleles to the child, i.e., the two parental alleles have equal probability of transmission, implying that the two ordered genotypes ( A i , A j ) and ( A j , A i ) both have probability p i p j when i ≠ j ; and 4) There is random mating between parents. This assumption finally implies that thus describing the (ordered) genotype distribution in the parental population. Note that when written as ordered genotypes, we can write P ( M , F , C ) = P ( M , F ) by notational convention, so that ( 1 ) actually describes the full expected distribution of triad genotypes in an unselected group. When analyzing ART triads, the validity of these assumptions can be checked using the non-ART triads. Below, we show below how deviations from these assumptions can be interpreted in terms of risk of ART use. This formulation covers any multi-allelic locus with H different alleles. Formulating the models for multi-allelic loci can be convenient, in particular to handle haplotypes constructed from a sequence of neighboring SNPs [ 29 ]. However, here it suffices to consider diallelic SNPs, with alleles A 1 and A 2 and where p 1  +  p 2  = 1. The A 1 allele will be designated as the “reference allele” and A 2 as the “effect allele”. In situations where we aim to estimate a relative risk parameter, we typically set RR 1  = 1 for the A 1 allele, and let RR 2 be the relative risk associated with the A 2 allele, a parameter to be estimated. For instance, at a C>T-SNP one might chose A 1  =  C as the reference allele and set A 2  =  T as the effect allele. In the current analyses, we typically define the major allele as the reference allele and the minor allele as the effect allele. We let x m and x f denote the parental allele dosages, so that x m ∈ { 0 , 1 , 2 } counts the number of A 2 alleles carried by the mother, and similarly with x f for the father and x c for the child. To model and estimate the parental genotype distribution P ( M , F ), we can write it using a log-linear formulation. A log-linear model assumes independent Poisson distributed counts of triads of each parental genotype combination, using a log-linear formulation for the expected numbers of each combination. As long as an intercept is included in the model, this is equivalent to writing it as a multinomial distribution, where the log probabilities follow a linear model: Here, α 0 can be regarded as a scaling constant that ensures the probabilities sum up to 1. Also, α = log ( p 2 / p 1 ) relates ( 2 ) to the multiplicative formulation in ( 1 ). Of particular note is that the linear dosage term α x m signifies that maternal alleles are in HWE, and similarly for paternal alleles, although this restriction can easily be lifted. Furthermore, the common regression coefficient α for both x m and x f means that mothers and fathers derive from populations with the same allele frequencies, and the fact that there is no interaction term between the mother ( x m ) and the father ( x f ) equates to assuming random mating between the two. If a given locus is involved in determining the risk of ART use for a couple, the parental genotype distribution P ( M , F | A ) among ART couples will be shifted relative to the unconditional distribution. More specifically, using Bayes’ formula, we have where the prevalence of ART, P ( A ), can be considered a normalizing constant. We thus see that the probability of ART use, P ( A | M , F ) , determines how the conditional distribution P ( M , F | A ) is shifted relative to P ( M , F ). In our previous work [ 24 ], we modeled the probability of ART use, P ( A | M , F ) , using a logistic regression model based on the genotypes of the individual parents, and combinations, i.e., interactions, of these. In the current setting, we use a log-linear model for the risk of ART, parameterized as where x m and x f are the dosage variables defined above, and β 0 is a scaling constant that makes the probabilities sum to 1. Note that the interaction effect β m f can be interpreted as an “incompatibility” between parental alleles, in that it impacts the probability of ART use above and beyond the individual contributions from each of the parents, as measured by β m and β f , respectively. The logistic model [ 24 ] used the non-ART triads as a baseline and compared the risk of ART use across different genotypes. However, in the current model, we can estimate parameters, in particular the interaction parameter β m f , only based on the ART triads. This strategy is similar to an established approach, case-only analyses, to estimate gene-environment (GxE) interactions. In case-only analyses, under the assumption of independence of genes and environment in the population, GxE effects can be estimated [ 35 , 36 ]. In our setting we can achieve this by applying ( 2 ) and ( 4 ) to ( 3 ), thus deriving the log-linear formulation with β 00 = α 0 + β 0 − log ( P ( A ) ) . It should be noted that when we are looking at ART families alone, using ( 5 ), we cannot estimate β m and β f directly; only the difference between them can be estimated by calculating ( α + β m ) − ( α + β f ) = β m − β f . In other words, when studying ART triads, we cannot infer whether the increased risk of ART use is attributable the maternal genotype, the paternal genotype, or both. It is also clear that estimating the difference β m − β f hinges on the assumption of a common α parameter for x m and x f from ( 2 ). Akin to a case-only study, the interaction effect β m f can be estimated based on the ART triads alone. However, note that this interaction effect may be aliased with non-random mating between parents. That is, if there is a correlation between parental genotypes, it would show up as a interaction term in ( 2 ). The above model description essentially covers Steps 1 and 2 in Figs 1 and 2 . Step 3, however, involves fetal genes, since alleles associated with early fetal survival would exhibit unexpectedly high frequencies in surviving children. Assuming that fetal survival during pregnancy depends both on fetal and maternal genes, the standard Haplin model [ 29 ] can be described as follows. Let S denote the outcome that the fetus survives pregnancy, i.e., we now consider S as the fetal “case” definition in a triad analysis. We then have for case-parent triads Here, P ( S ) is again just a normalizing constant. The standard multiplicative model formulation for the penetrance used in Haplin is then where B is a baseline risk, R R j , R R l are fetal relative risks associated with the fetal A j , A l alleles, respectively, and R R i ( M ) , R R j ( M ) are relative risks associated with the maternal alleles A i , A j , respectively. Since A 1 is assigned the role of reference allele, only RR 2 and R R 2 ( M ) need to be estimated; RR 1 and R R 1 ( M ) are both set to 1. More details about model implementation and interpretation, and the application of log-linear models to case-parent triads, are available in previous publications [ 28 , 29 , 37 ]. Note that if none of the alleles confers any risk, then P ( S | M , F , C ) is constant, and so that the parental genotype distribution would be the same among couples that experience a successful birth as among all couples. Transforming ( 7 ) to log scale, we obtain for suitable γ -parameters, and applying ( 2 ) and ( 8 ) to ( 6 ) we obtain with γ 00 = α 0 + γ 0 − log ( P ( S ) ) . Again, only the difference ( α + γ m ) − α = γ m can be estimated. We now combine the two model approaches into a joint model for ART parents with a surviving fetus. We again focus primarily on modeling the genotype distribution of the ART triads, this time conditioning on both A and S , i.e., P ( M , F , C | S , A ) . As before, we can express this as Note that within ART triads with a surviving fetus, P ( S , A ) is again just a normalizing constant in this expression. The penetrance model for fetal survival, P ( S | M , F , C , A ) , can be approached as in ( 8 ) above, except with the understanding that the risk parameters γ c and γ m are estimated conditional on ART use, i.e., it measures genetic effects on fetal survival in triads where ART is being used. The P ( A | M , F , C ) and P ( M , F ) parts can be modeled as described above. By combining ( 2 ), ( 4 ), and ( 8 ), we obtain the following log-linear model with δ = α 0 + β 0 + γ 0 − log ( P ( S , A ) ) . To estimate model parameters, we again assume the observed count of each type of ( M , F , C ) triad follows independent Poisson distributed variables, with expected values proportional to P ( M , F , C | S , A ) in ( 9 ). As discussed above, to ensure identifiability of parameters, some parameters in ( 9 ) have to be restricted. While the difference ( α + β m + γ m ) − ( α + β f ) = β m + γ m − β f can be estimated when again conditioning on the common parameter α , the parameters β m , γ m , and β f cannot be distinguished. However, the two most important parameters in our analysis, γ c and β m f , can be identified directly. These correspond to the two study outcomes described previously: the parental interaction parameter β m f captures the effect on the parental risk of requiring ART to conceive, while γ c reflects the effect of the fetal genotype on fetal survival. The completed model ( 9 ) thus incorporates all three steps in Fig 1 , and under suitable assumptions and parameter restrictions, the effects of the three steps can be separated. Violations of model assumptions may influence the effect estimates to varying degrees. The fetal effect is estimated conditional on parental genotypes and is therefore only modestly affected by population stratification [ 28 ]. In contrast, the parental interaction effect is more vulnerable to bias, particularly from population stratification and assortative mating. If two individuals in a couple are more genetically similar than expected by chance, this similarity may be misattributed by the model as a parental interaction effect. Furthermore, in relation to Step 3, it should be noted that while the model includes the fetal survival effect in ART triads, as expressed in ( 8 ), it does not explicitly model the indirect selection effect on parental alleles caused by possible previous pregnancy losses for the same parents; the C in the formulas denotes the genotype of the successful ART birth, not any previously lost (and thus unobserved) siblings. It should also be recognized that although the model in ( 9 ) is derived within the context of only ART triads, comparable analyses can also be conducted within non-ART triads following similar reasoning. Note that in all the above model descriptions, the dosage variables were added in a linear fashion, i.e., they are assumed to have a multiplicative dose-response effect. However, in many situations, it would make more sense to allow a free response model, i.e., a model that allows different relative risks for a single- and a double-dose of the effect allele. To achieve this, we can replace the dose variables by their components, for instance by replacing x c with x c , 1 + x c , 2 . Here, x c ,1 is 0 if the first allele of the child is A 1 , and 1 if it is A 2 , and x c ,2 similarly defined for the second allele of the child. Then, including a within-genotype allelic interaction in ( 8 ), the model γ c x c , 1 + γ c x c , 2 + γ c , 12 x c , 1 x c , 2 estimates the parameter γ c , 12 as a deviation from expectations of the multiplicative model; it thus allows for, for instance, recessive and dominant models to be fitted. Another example is if both x m and x f are split into two parts in (2) ; the corresponding model with within-genotype interactions allows deviations from the HWE assumption. In line with our previous modeling of parental interaction effects [ 24 ], we implemented four different interaction models here: multiplicative, dominant, complementary, and threshold-based. The multiplicative interaction effect model was implemented as described in ( 9 ), based on the interaction model in ( 4 ). The dominant parental interaction effect model requires both parents to carry at least one copy of the minor allele, but the effect of one vs. two copies of the minor allele is assumed to be equal. To illustrate, we then define x ~ m as 0,1,1 for the maternal genotypes A 1 A 1 , A 1 A 2 , A 2 A 2 , respectively, and similarly for x ~ f . Thus, the interaction term x ~ m x ~ f equals one in the case where both parents carry at least one copy of A 2 , otherwise it is zero. These definitions of x ~ m and x ~ f would then be used in the interaction model ( 4 ) instead of x m and x f , and ( 4 ) would again be inserted in ( 9 ), keeping x m and x f separate from the new x ~ m and x ~ f . With a complementary parental interaction effect, we are targeting parents who are homozygous for the opposite alleles, i.e., at a specific locus, one parent is homozygous for the minor allele while the other is homozygous for the major allele. Finally, the threshold-based parental interaction effect assumes that the total number of minor alleles carried by both parents have an added effect on the outcome if they exceed a given threshold, in this setting equal to two. Similarly to the dominant model, both the complementary and the threshold-based models are implemented by modifying the interaction model in ( 4 ) accordingly before inserting it into ( 9 ). For more details on the parental interaction models and the specific coding of the complementary and threshold-based models, see Skodvin et al. [ 24 ]. Although not immediately evident from ( 9 ), the models assuming log-linear dose relationships, such as those containing the terms γ c x c and β m f x m x f using the dosage variables, exhibit minimal aliasing between parental and fetal components, as explained in detail previously [ 28 ]. In our setting, fetal and parental effect estimates start impacting one another particularly for the complementary and threshold parental models in combination with a free model for the fetal genes. Consequently, these are the combinations we have focused on in our current analyses. The complementary and threshold models have some overlap but also differ in ways that are important to consider when interpreting the parental interaction effect estimates. Parents with the complementary genotype combination will also meet the threshold of two minor alleles as required by the threshold model. However, such couples likely represent only a minority among those identified by the threshold model. While the complementary model favors genetically dissimilar parents, the threshold model assigns larger interaction effects as the total number of minor alleles carried by both parents increases, regardless of whether these alleles come from similar or dissimilar genotypes. We implemented the log-linear model described in ( 9 ) as an extension of the log-linear model already available in Haplin. By using maximum likelihood estimation and applying the expectation-maximization (EM) algorithm, missing data can be imputed, which means that Haplin also allows for analyzing incomplete triads, such as mother–child or father–child dyads. It also allows probability weighting of all possible configurations in triads with ambiguous transmission of alleles from parents to child, as well as imputation and weighting of possible haplotype configurations across multiple SNPs. For a detailed description of Haplin, see the original article [ 29 ]. We used data from MoBa, a nationwide pregnancy cohort to which pregnant women were recruited from 1999 to 2008. Expectant mothers were invited to join the study at the time of their routine ultrasound examination, typically around week 17 of gestation, and the participation rate was 41% [ 38 ]. From 2001, fathers were also invited to the study. Currently, there are approximately 95,200 mothers, 75,200 fathers, and 114,500 children in the cohort. Blood samples were collected from the mothers and fathers at the time of the routine ultrasound examination and from the mothers and their newborn’s umbilical cord at birth [ 39 ]. Due to the recruitment strategy of MoBa, the data are by default restricted to parents who were able to achieve a successful fertilization and carry the pregnancy to at least approximately week 17 of gestation. Hence, the use of ART serves as a key indicator of fertility challenges. Information on ART use was retrieved through data linkage with the Medical Birth Registry of Norway (MBRN) [ 40 ]. ART use was classified based on records indicating whether the pregnancy was achieved through in vitro fertilization or intracytoplasmic sperm injection. We included siblings in the analyses, but excluded multiple pregnancies. The two study samples eligible for analysis were as follows: for the non-ART sample, we had 42,047 mother-father-child triads; for the ART sample, we had 842 mother-father-child triads, 411 mother-child dyads, and 83 father-child dyads. Our data included a few instances of late stillbirths, with 18 occurring among non-ART pregnancies and fewer than 5 among ART pregnancies (exact number not provided to ensure privacy). All were reported to have occurred after gestational week 22. Fig 1 shows that the group of ART parents comprises four subgroups from Step 1 based on parental fertility status: infertile mother, infertile father, both parents infertile, or neither parents infertile. Although complementary information on infertility-related conditions is available in MBRN for some individuals, our data do not provide complete information to fully classify these subgroups. Nevertheless, our modeling strategy accounted for the main maternal and paternal effects, as demonstrated in ( 4 ). In other words, these effects are not ignored, even though they cannot be directly estimated within an ART triad design, as noted previously. Importantly, the maternal–paternal interaction effect—one primary focus of this study—can be estimated. An extensive effort has been invested in genotyping as many individuals as possible in the entire MoBa cohort, with the process undertaken in multiple batches over the years [ 41 ]. Currently, the genetic data comprise approximately 77,000 mothers, 53,000 fathers, and 76,000 children. The genetic data have undergone thorough quality control as detailed by Corfield et al. [ 41 ]. We further processed the genetic data using PLINK v1.90 [ 42 ], applying a minor allele frequency (MAF) threshold of less than 1% and pruning SNPs in strong linkage disequilibrium ( r 2 >0.95), leaving approximately 2,300,000 SNPs for the current analysis. We used Haplin [ 30 , 31 ], implemented in R v4.0.4 [ 43 ], to analyze fetal effects on fetal survival, initially without adjusting for parental interaction effects. We then reanalyzed the fetal effects using the extension of the Haplin models including adjustments for parental interaction effects. The fetal effects were estimated by separate parameterizations of the single- and double-dose effects of the effect allele, with complementary and threshold-based parameterizations of the parental interaction effect. The analyses were conducted separately in the non-ART and ART samples. The effects were quantified as relative risk (RR) estimates, with 95% confidence intervals (CIs) and p -values. Further, we calculated the false discovery rate (FDR) q -values by applying the Benjamini-Hochberg method [ 44 ]. The Haplin framework also includes a test for Hardy-Weinberg equilibrium (HWE). All visualizations of the results were created using the R packages Haplin [ 30 ] and ggplot2 [ 45 ].

(XLSX) (XLSX) (XLSX) Panel A displays the negative log 10 -transformed p -values, and the red line indicates the Bonferroni-corrected genome-wide significance threshold of 5 × 10 −8 . Panel B shows the corresponding quantile-quantile (QQ) plot, along with the genomic inflation factor ( λ ). (TIF) Panels A and C display the negative log 10 -transformed p -values, and the red line indicates the Bonferroni-corrected genome-wide significance threshold of 5 × 10 − 8 . Panels B and D show the corresponding quantile-quantile (QQ) plots, along with the genomic inflation factor ( λ ). (TIF) Panels A and C display the negative log 10 -transformed p -values, and the red line indicates the Bonferroni-corrected genome-wide significance threshold of 5 × 10 − 8 . Panels B and D show the corresponding quantile-quantile (QQ) plots, along with the genomic inflation factor ( λ ). (TIF) Panels A and C display the negative log 10 -transformed p -values, and the red line indicates the Bonferroni-corrected genome-wide significance threshold of 5 × 10 − 8 . Panels B and D show the corresponding quantile-quantile (QQ) plots, along with the genomic inflation factor ( λ ). (TIF) Panels A and C display the negative log 10 -transformed p -values, and the red line indicates the Bonferroni-corrected genome-wide significance threshold of 5 × 10 − 8 . Panels B and D show the corresponding quantile-quantile (QQ) plots, along with the genomic inflation factor ( λ ). (TIF) Panels A and C display the negative log 10 -transformed p -values, and the red line indicates the Bonferroni-corrected genome-wide significance threshold of 5 × 10 − 8 . Panels B and D show the corresponding quantile-quantile (QQ) plots, along with the genomic inflation factor ( λ ). (TIF) The green arrows indicate the direction of transcription for the genes in the region, while recombination rates are represented by blue lines. The color gradient reflects the degree of linkage disequilibrium, with associated r 2 values. (EPS) The parental interaction effect are shown for the complementary model (panel A ) and the threshold-based model (panel B ) in the ART sample (blue) and the non-ART sample (red). The negative log 10 -transformed p -values are displayed, and the red lines indicate the Bonferroni-corrected genome-wide significance threshold of 5 × 10 − 8 . (TIF) Panels A and C display the parental complementary interaction effects, while panels B and D display the parental threshold-based interaction effects. The shaded area represents a 95% confidence interval band for the null hypothesis p -value distribution. SNPs with MAF lower than 5% are excluded. (TIF) Statistical power was estimated as the proportion of analyses in which the parental complementary interaction effect was significant at the Bonferroni-corrected threshold of 5 × 10 − 8 . (EPS) (EPS)