Polygenic scores for dizygotic twinning: insights into the genetic architecture of female fertility

TL;DR

This study explores the genetic factors behind dizygotic twinning and how they relate to female fertility and infertility.

Contribution

The study introduces a polygenic score for dizygotic twinning and links it to fertility traits and treatment outcomes.

Findings

Dizygotic twinning has low polygenicity, with only 0.20% of SNPs having a nonzero effect.

The DZT polygenic score explains 1.6% of variance in DZT liability and distinguishes between natural and medically assisted pregnancies.

DZT shows strong genetic correlations with anovulatory infertility and PCOS, suggesting a genetic spectrum of fertility.

Abstract

Natural dizygotic twinning (DZT) results from hyper-ovulation and is an indicator of female fertility. However, some traits linked to DZ twinning are also associated with infertility. We examined the relationship between DZT and female (in)fertility using recent GWAS findings.

We investigated the genetic architecture of DZT and compared polygenic scores (PGS) for DZT between mothers of naturally conceived DZ twin pregnancies and mothers who required fertility treatments (MAR) in the Netherlands Twin Register (NTR) and the Norwegian Mother, Father, and Child Cohort Study (MoBa). We also calculated genetic correlations between DZT and seven fertility related traits.

DZT has a low polygenicity, with only 0.20% of SNPs estimated to have a nonzero effect. The DZT PGS explains 1.6% of variance in DZT liability, and we observe an odds ratio of 2.29 between the first and the tenth PGS deciles. The DZT PGS distinguishes between mothers of naturally conceived pregnancies and mothers who received MAR and is associated with a shorter time to pregnancy in mothers of singletons. The lowest PGSs were observed for mothers who received hormonal ovulation induction, indicating maternal fertility issues. DZT showed genetic correlations with anovulatory infertility (rg =  − 0.698) and PCOS (rg =  − 0.278), and endometriosis (rg = 0.279).

Female fertility appears to exists on a genetic spectrum, with anovulation/infertility at one end and DZT at the other. Results suggest that the DZT PGS can be of added value to evaluate female fertility and be incorporated in clinical practice in the future.

The online version contains supplementary material available at 10.1007/s10815-025-03515-4.

Figures

The variance explained on the liability scale by the polygenic score (PGS) for dizygotic (DZT) twinning in 1635 cases and 2797 controls for each fraction of top SNPs from the 2024 DZT GWAS included in the PGS

Odds ratios for dizygotic twinning (DZT) in the 7178 Netherlands Twin Register (NTR) participants by PGS decile. Odds ratios and 95% confidence intervals were estimated with logistic regression for both A being of mother of a naturally conceived DZ twin and B being a mother of a DZ twin who conceived after medically assisted reproduction (MAR) treatments. The reported odds ratios are relative to the fifth decile. The points represent odds ratios, and the bars represent the lower and upper 95% CI of the odds ratios

The mean PGS for dizygotic (DZ) twinning in mothers of naturally conceived DZ, controls, and monozygotic (MZ) twins from the A and C Netherlands twin register (NTR) and B and D the Norwegian Mother, Children, and Father (MoBa) cohort studyTable 2The comparison of the polygenic score for dizygotic (DZ) twinning in mothers of DZ twins with mothers of MZ twins and singletons from the Norwegian Mother, Children, and Father (MoBa). The controls include the mothers of naturally conceived singletons. All the results are shown with the mothers of naturally conceived DZ twins as the reference group. MA

The mean PGS for dizygotic (DZ) twinning in mothers of naturally conceived DZ twins and mothers from DZ twins that received medically assisted reproduction (MAR) treatments from A the Netherlands Twin Register (NTR) and B the Norwegian Mother, Children, and Father Cohort Study (MoBa)Fig. 5The mean PGS of the DZ medically assisted reproduction (MAR) group of the NTR is further split between ovulation induction (OI) and IVF and ICSI

Full text

Introduction

Giving birth to naturally conceived dizygotic (DZ) twins starts with poly-ovulation, followed by multiple fertilization and successful implantation of at least two embryos [1]. DZ twinning (DZT) occurs more often in some families than in others and is partially genetically influenced through the mothers’ side, which is what we would expect for a trait that is dependent on the mother’s physiology [1–4]. Two genome-wide association studies (GWAS) report multiple independent loci associated with DZT [5, 6]. The SNP-based heritability was estimated to be 2.4% on the liability scale in the most recent GWAS meta-analysis of DZT [6]. The heritability estimate for DZ twinning is 28%, based on a reanalysis of sister-pair data [4]. An analysis of seven large multi-generation pedigrees from West Africa, fin-de-siècle French Jewish populations, Canadians, and the French royal family obtained an estimate of 8 to 20%, without distinguishing between monozygotic and dizygotic twins [4, 7].

DZT presents as an “atypical” complex trait in genetic studies. The phenotype definition is critically important as illustrated by the first GWAS for DZT, when a genetic signal was observed only after removing a small subset (N = 241) of mothers of DZ twins who had conceived twins after medically assisted reproduction (MAR) [8]. Such a strong effect of phenotype misclassification hints at opposing genetic effects associated with the natural conception of DZ twins versus mothers with subfertility or infertility. Also, in the most recent meta-analysis, the genomic variants that were identified generally are located in the coding regions of the genome, as opposed to non-coding regions that are often found in GWAS [9, 10]. Predicted causal genes are involved in endocrine processes of female fertility including follicle stimulating hormone subunit B (FSHB), gonadotropin releasing hormone 1 (GNRH1), and follicle stimulating hormone receptor (FSHR) [6, 11].

It has been proposed that DZT can be viewed as an indicator of high female fertility and that research into DZT can thus aid in understanding female fertility and infertility in general [12] Paradoxically, several characteristics associated with DZT are associated with infertility [13–15]. Infertility is defined as being unable to conceive after 12 consecutive months of trying and has been estimated to affect one in six couples worldwide [16]. Several underlying conditions may lead to female infertility including polycystic ovary syndrome (PCOS) and endometriosis, but for at least 15% of couples, the underlying reason for infertility remains unclear [17]. MAR treatments for infertility have been available since the late 1960 s and can result in increased dizygotic twinning rates for two reasons: [1] increase in hyper ovulation as a result of hormonal ovulation induction (OI) mimicking the natural DZT mechanism [1, 18] or [2] the transfer of multiple fertilized oocytes after in vitro fertilization (IVF) cycles to the uterus, although the transfer of more than one fertilized oocyte has become less common in recent years [19, 20].

Our objectives were to investigate the genetic architecture of DZT and to test the hypothesis that DZT can be regarded as an indicator of female fertility. We first investigated the polygenicity of DZT by employing multiple tools [21–24]. Next, we compared the polygenic score (PGS) for DZT between twin and singleton mothers who conceived naturally and mothers who received MAR treatments. We tested our hypothesis in the Netherlands Twin Register (NTR) and in the Norwegian Mother, Father and Child Cohort Study (MoBa). We extended the PGS analyses by assessing whether the DZT PGS was associated with the time to achieve pregnancy in singleton mothers. To further understand the relation between (sub)fertility and DZT, we estimated the genetic correlations between DZT and seven fertility phenotypes [21, 22].

Methods

Participants and phenotypes

NTR

The Netherlands twin register (NTR) is a national twin-family register established in the late 1980 s [25, 26]. A total of 4381 genotyped mothers were included and categorized into four groups: mothers of natural DZ twins (N = 1635), mothers of natural MZ twins (N = 2107), MAR MZ mothers (N = 125), and MAR DZ mothers (N = 514). MAR was defined as the use of IVF, intracytoplasmic sperm injection (ICSI), and OI. The DZ MAR mothers were further stratified into mothers of DZ twins who received hormonal induction ovulation without IVF (OI; N = 252) and the IVF/ICSI DZ mothers (N = 262) who underwent IVF or ICSI. Mothers were excluded from the study if they were not the biological mother (N = 92), if they were a mother of both MZ and DZ twins (N = 65), or if they were mothers of MZ twins and received OI treatments (N = 27). We included 2797 fathers of twins as controls. Given that the fathers were screened on not being a multiple themselves and that there is no evidence for assortative mating for twinning, no confounding effects were expected by using fathers of twins as controls [27]. Mothers without information on MAR were classified as non-MAR when the twins were born before 1970, as IVF was introduced in the Netherlands in the late 1970 s and clomid (a hormonal substitute used to induce ovulation) in the late 1960 s [28]. Surveys provided information on maternal age at birth of twins, age first birth, height, BMI before pregnancy, number of older children, and smoking before the twin pregnancy [15]. We determined time to pregnancy (TTP) for 2298 mothers of naturally conceived twins. TTP was determined in four categories: 0–2 months, 3–5 months, 6–12 months, and more than 12 months.

MoBa

The Norwegian Mother, Father and Child Cohort Study (MoBa) is a population-based pregnancy cohort study conducted by the Norwegian Institute of Public Health [29]. Participants were recruited from maternity wards all over Norway from 1999 to 2008. The women consented to participation in 41% of the pregnancies. Blood samples were obtained from both parents during pregnancy and from mothers and children (umbilical cord) at birth [30]. The cohort includes approximately 114,500 children, 95,200 mothers and 75,200 fathers. The current study is based on version 12 of the quality-assured data files released for research in January 2019.

Data on the participating mothers were linked with data from the Medical Birth Registry of Norway (MBRN; a national health registry containing information about all births in Norway), to assess their reproductive history. We included 1299 mothers of DZ twins and 64,810 mothers of singletons. Mothers were excluded from the singleton group if they had an earlier or later twin pregnancy (N = 2110). We categorized mothers of twins into three groups: mothers of natural DZ twins (N = 658), mothers of MAR DZ twins (N = 328), and mothers of natural MZ twins (N = 313). Mothers of singletons were categorized as mothers of natural singletons (N = 63,422) or mothers of MAR singletons (N = 1388). The mothers of naturally conceived singletons were the control group in the MoBa analyses. Mode of conception was reported by record linkage to the MBRN, which contains information on IVF/ICSI use as reported from the fertility clinics, or by self-report in MoBa. Other information on maternal age, reproductive history and maternal body composition and TTP was obtained from MoBa questionnaires or from the MBRN [15, 31]. We determined TTP for 42,285 mothers of naturally conceived singletons and TTP ranged between 0 (unplanned pregnancy) and 12 or more months.

PGS calculation

NTR

The complete genotyping and imputation information from the NTR can be found in Data S1. The PGSs in the NTR sample were based on the discovery GWAS meta-analysis for DZT of Mbarek et al. [6] A leave-one-out (LOO) meta-analysis was performed, removing all NTR summary statistics from the meta-analysis and rerunning the meta-analysis as detailed in Mbarek et al. (2024). In short, the NTR GWAS results were removed and the GWAS meta-analysis for DZT was run for the remaining cohorts using the METAL software that performs a fixed effect meta-analyses [32].We retained variants for which the effect allele frequency (EAF) was 0.01 ≤ EAF ≤ 0.99. Variant EAF and effect sizes of the 1000 genome variants of the PRS summary statistics were allele and strand aligned with the imputed NTR genetic data. Discovery variants that could not be aligned due to allele mismatches before or after strand flipping were discarded. The processed LOO summary statistics were taken as input for the LDpred 0.9 softwareV1. For estimating the target LD structure, we took a selection of 2500 unrelated individuals in the NTR database and selected a set of well-imputed variants in the NTR sample [21, 33].

The parameter ld_radius was set by dividing the number of variants in common (from the output of the coordination step) by 12,000, ensuring an LD window of 500 kb (on average) on either side of each SNP [34]. For the coordination step, we provided the median sample size as input value for N. For the LDpred step, we applied the following thresholds for fraction of variants with a nonzero effect (in addition to the default infinitesimal model): –PS = 0.5, 0.3, 0.2, 0.1, 0.05, 0.01, 0.005, 0.003, and 0.001.

In addition, we calculated a PGS in NTR based on the 26 unique significant SNPs that were also present in the Human Genome Diversity Project and the 1000 Genomes Project, from the most recent DZT GWAS (Table 2 from Mbarek et al., 2024). SNPs were selected from the LOO summary statistics, and the PGS were calculated using the –score function of PLINK version 1.90.

MoBa

The complete genotyping and imputation information from MoBa can be found in Data S1. In MoBa, we calculated the PGS in a similar way as in NTR. The full GWAS summary statistics, including the NTR, were taken as input for the LDpred 0.9 software. For estimating the target LD structure, the same sample was used [21, 33]. The parameter ld_radius was set by dividing the number of variants in common (from the output of the coordination step) by 12,000. Again, for the coordination step, we provided the median sample size as input value for N. For the LDpred step, we applied the 0.001 threshold for fraction of variants with nonzero effect. We chose to only include the 0.001 threshold for MoBa, based on the variance explained results from the NTR PGSs.

Statistical methods

Polygenicity estimation

We employed the GENESIS tool as described by Zhang et al. (2018) for estimating the genetic effect-size distribution based on the summary statistics of the 2024 GWAMA of DZT [6, 23]. We first tested whether a two-components model (M2) or a three-components model (M3) fitted our data best as described in Zhang et al. (2018). The M2 model entails a model where the effects of the nonzero effect SNPs will follow a normal distribution, while M3 implies that a mixture of SNPs with a small effect and a larger effect will influence the trait. We based the decision of the best fitting model on the BIC (Bayesian Information Criterion) and the concurrence between the observed and expected lambda as indicated in the QQplot.

We also estimated the proportion of SNPs with nonzero effects as an indication of the polygenicity in the SBayesS tool [24]. For the SbayesS analyses, we used the standard parameters and a LD matrix based on UK Biobank (UKB) provided on the SbayesS webpage (https://cnsgenomics.com/software/gctb/#Download).

Variance explained by PGS on the liability scale

All analyses were executed in R4.3.2. For each set of weighted effect sizes (infinitesimal model and additional threshold, see above) we calculated the NTR PGSs over all genotype data sets. To select the optimal threshold value of the fraction of variants with nonzero effect, we calculated the variance in DZT explained on the liability scale in 1635 mothers of naturally conceived DZ twins and 2797 controls employing the methods described by Lee et al. [35].

Odds ratios

For all 7178 participants of the NTR, we determined in which decile their PGS was located. We then performed a logistic regression (glm function) for being a mother of a naturally conceived DZT and separately for being a MAR DZT mother on the PGS deciles. From this regression, we determined the odds ratios and corresponding 95% confidence interval for each of the deciles using the exp function in rbase.

PGSs comparisons in NTR and MoBa

We performed logistic regression with PGS as predictor and case–control as outcome variable, corrected for relatedness by adding the family structure as a cluster (function geeglm) as both the NTR and MoBa sample included sister pairs. The family structure was defined by sisters sharing the same family ID. To correct for population stratification and genotyping/batch effects, the first 10 PCs were added as covariates as well as two dummy-coded variables for three genotyping platforms in NTR. In MoBa, genotyping occurred in 24 batches, which consist of varying genotyping platform and facility (Supplement Data 1). The batch and genotyping platforms were therefore included as a covariate in the MoBa analysis as a factor. To investigate the effect of DZT related traits, we performed a second round of analyses that included covariates that were previously associated with DZT, namely, maternal age at birth, height, BMI, smoking, and reproductive history [15, 36].

We performed the following series of analyses in NTR and MoBa. Bonferroni adjusted significance thresholds were 7.1 × 10^−3^ (0.05/7) in NTR and 0.01 (0.05/5) in MoBa.

1. Mothers of naturally conceived DZ twins compared to controls. In NTR, controls are fathers of twins (non-twins themselves), and in MoBa, controls were mothers of naturally conceived singletons.
2. Mothers of naturally conceived DZ twins compared to mothers of naturally conceived MZ twins.
3. Mothers of naturally conceived DZ twins compared to mothers of MAR DZ twins. We further stratified the MAR DZ mothers in mothers who received OI treatment and IVF/ICSI mothers in the NTR.
4. Controls compared to mothers of naturally conceived MZ twins and mothers of MAR MZ twins.
5. In MoBa only: Mothers of naturally conceived singletons compared to mothers of MAR singletons.
6. Regression analyses with the PGS in mothers of naturally conceived twins (NTR) and singletons mothers (MoBa), and the TTP. TTP was analyzed both as a continuous variable and as a categorical variable stratified per months.

Genetic correlations

We applied LD Score Regression (LDSC; version v1.0.1; [22, 37]) to calculate the genetic correlation of DZT with number of children ever born (NEB) and childlessness (CL) [38], PCOS [39], anovulatory infertility [40], endometriosis [41], miscarriages [42], and male infertility [43]. The summary statistics of endometriosis were downloaded from the GWAS catalog, while the other summary statistics were made available by the authors of these studies. Before the genetic correlations were calculated, the SNPs for each of the summary statistics were merged with the w_hm3.snplist to ensure that the alleles listed in the summary statistics files match, which was downloaded from the LDSC wiki (https://github.com/bulik/ldsc/wiki/).

Results

Polygenicity of DZ twinning

To investigate the polygenicity of DZT, we used two genetically informed tools. We employed GENESIS (https://github.com/yandorazhang/GENESIS) for predicting the genetic effect-size distribution based on the summary statistics of the 2024 GWAMA of DZT [6, 23]. We fitted both the M2 and M3 models that GENESIS provides. The M2 models entail models where a large number of SNPs will have a small effect on the trait, while M3 implies that a mixture of SNPs with a small effect and a larger effect will influence the trait. By assessing BIC and the fit of the QQ plots, Figure S1 shows that the M3 model was best suited implying that DZT is influenced by a mixture of SNPs with a small effect and a larger effect a mixture of SNPs with a small effect and a larger effect. We calculated an effective sample size of discovery GWAS for DZT by the described formula: (8265 + (26,252/2)) × 682,000/716,517) = 20,361. Mothers of DZ twins are (N = 8265) counted as one case, while DZ twins (N = 26,252) are counted as half cases given that the twins carry half of their mothers genes. The total predicted number of susceptible SNPs for DZT is around 2244 (SD = 1498), and the predicted maximum SNP heritability based on these 2244 SNPs was 7.2% (SD = 2.4%). Of the 7.2% SNP heritability, 17.2% (SNP heritability = 1.2%; SD = 0.6%) is predicted to be explained by a cluster of only 17 SNPs with the largest effect size. The results of GENESIS also show that a sample size of approximately 400,000 participants (200,000 cases and 200,000 controls; Figure S2) would be required to identify 80% of the total SNP-based heritability.

Next, we ran SbayesS to evaluate the polygenicity of DZT. SbayesS estimated the proportion of SNPs with nonzero effects as an indication of the polygenicity [24] https://cnsgenomics.com/software/gctb/). We used a LD matrix based on UK Biobank (UKB) provided on the SbayesS webpage (https://cnsgenomics.com/software/gctb/#Download) and the summary statistics of the DZT GWAS by Mbarek et al. (2024) [6]. SbayesS estimated that about 0.198% (SE = 4.00x^−5^) or approximately 5483 (SD = 939) of SNPs will have a nonzero effect with regard to DZT. The two estimates of the susceptible SNPs from GENESIS and SbayesS are not significantly different from each other based on the surrounding confidence intervals, indicating robustness in the estimates. SbayesS estimates the SNP based heritability of DZT based on the 2024 GWAS at 0.8% (SE = 9.19x^−5^).

Variance explained by the DZT PGS

To investigate the optimal fraction of SNPs to include in the PGS, we calculated the explained variance in NTR on the liability scale in a sample of 1635 natural mothers of DZ twins and 2797 controls. We first calculated PGS using the LDpred 0.9 softwareV1 (Methods) with several fraction of SNPs (f) with nonzero effects (infinitesimal, 0.5, 0.3, 0.2, 0.1, 0.05, 0.01, 0.005, 0.003, and 0.001). Figure 1 shows that explained variance range of 0.71–1.60% on the liability scale (Table S1). The highest variance explained was observed for f = 0.001, and thus, this was the fraction selected to calculate the PGS in NTR and MoBa. Figure 1 also includes a PGS solely based on 26 genome-wide significant SNPs which in the 2024 GWAS were included in a PGS to explore global differences in DZT. These SNPs were selected on being a unique hit as indicated by FUMA and on being present in the Human Genome Diversity Project and the 1000 Genomes Project databases. Clearly, this PGS had the lowest explained variance (R^2^ = 0.71%), pointing to a need for more research in non-EU populations.Fig. 1. The variance explained on the liability scale by the polygenic score (PGS) for dizygotic (DZT) twinning in 1635 cases and 2797 controls for each fraction of top SNPs from the 2024 DZT GWAS included in the PGS

PGS analyses

Odds ratios

The PGS of all 4381 mothers of twins and 2797 controls in NTR were and split into 10 decile groups. We calculated odds ratios for giving birth to a naturally conceived DZ twin and giving birth to a DZ twin after MAR treatments in each decile compared to the first and the fifth deciles. Figure 2 A summarizes the odds ratio compared to the fifth decile, indicating a 1.67 increased risk of giving birth to a DZ twin in the tenth decile compared to the fifth decile (P = 2.91x^−5^). Between the first and tenth deciles, the increase in risk was 2.29 (P = 1.60x^−10^). Figure 2B indicates an 1.58 increased risk for needing MAR treatments to conceive a DZ twin in the first decile compared to the fifth decile (P = 0.024) and a 1.54 increased risk in the second decile compared to the fifth decile (P = 0.034). The odds of needing MAR treatments were not increased between the first and tenth decides.Fig. 2. Odds ratios for dizygotic twinning (DZT) in the 7178 Netherlands Twin Register (NTR) participants by PGS decile. Odds ratios and 95% confidence intervals were estimated with logistic regression for both A being of mother of a naturally conceived DZ twin and B being a mother of a DZ twin who conceived after medically assisted reproduction (MAR) treatments. The reported odds ratios are relative to the fifth decile. The points represent odds ratios, and the bars represent the lower and upper 95% CI of the odds ratios

Controls and mothers of naturally conceived DZ twins

We compared the DZT PGS between mothers of natural DZ twins and controls from the NTR using the PGS based on the f = 0.001 of the SNPs with a nonzero effect. Table 1 and S2 and Fig. 3A show that the mean PGS was significantly higher in the mothers of naturally conceived DZ twins compared to the controls (P = 6.94 × 10^−13^). In MoBa, the mothers of natural DZ twins also had a significantly higher mean PGS for DZT compared to controls (the mothers of naturally conceived singletons; P = 1.20 × 10^−4^; Table 2 and S3; Fig. 3B). Table 1. The comparison of the polygenic scores (PGS) for dizygotic (DZ) twinning in mothers of DZ with controls and other mothers of twins from the Netherlands Twin Register (NTR). The controls consist of fathers of twins. All the analyses are shown with the mothers of naturally conceived DZ twins as the reference group. MAR = medically assisted reproduction; OI = Ovulation inductionGroupN**P valueBeta (95% CI)Mean maternal age (sd)Mean age first birth (sd)Mean height (sd)Mean BMI before pregnancy (sd)Number of older children (sd)% Smoking everDZ Naturally1.635--30.00 (4.08)27.61 (4.02)170.16 (6.54)22.80 (3.24)0.92 (0.94)21.35%Controls2.7976.94 × 10^−13^ − 1.06 × 10^8^ (− 1.35 × 10^8^ to − 7.64 × 10^7^)31.89 (4.63)30.16 (4.59)182.33 (7.21)-^^1.04 (0.89)23.47%DZ MAR5141.33* × 10^−6^ − 1.16 × 10^8^ (− 1.63 × 10^8^ to − 6.79 × 10^7^)31.80 (3.75)30.76 (3.90)170.14 (6.22)23.13 (3.86)0.32 (0.52)11.91%DZ OI2522.00 × 10^−9^ − 1.92 × 10^8^ (− 1.28 × 10^8^ to − 2.56 × 10^8^)30.39 (4.07)29.25 (4.13)169.20 (6.04)23.02 (4.01)0.39 (0.56)13.89%DZ IVF/ICSI2620.0903 − 5.29 × 10^7^ (− 1.15 × 10^8^ to 9.58 × 10^6^)33.15 (3.44)32.21 (3.67)171.05 (6.39)23.24 (3.71)0.25 (0.49)9.92%MZ Naturally2.1072.09 × 10^−10^ − 9.94 × 10^7^ (− 1.31 × 10^8^ to − 6.81 × 10^7^)29.74 (4.19)27.89 (3.98)169.84 (6.34)22.98 (3.50)0.68 (0.84)16.50%MZ MAR1250.0884 − 7.25 × 10^7^ (− 1.58 × 10^8^ to 1.26 × 10^7^)34.18 (3.61)33.10 (3.49)171.20 (6.64)22.95 (2.46)0.37 (0.67)8.06%Italic inidcates a significant effect based on the alpha of 7.1 × 10-3^^Controls of the NTR only consist of males therefor BMI before pregnancy is not applicableFig. 3The mean PGS for dizygotic (DZ) twinning in mothers of naturally conceived DZ, controls, and monozygotic (MZ) twins from the A and C Netherlands twin register (NTR) and B and D the Norwegian Mother, Children, and Father (MoBa) cohort studyTable 2The comparison of the polygenic score for dizygotic (DZ) twinning in mothers of DZ twins with mothers of MZ twins and singletons from the Norwegian Mother, Children, and Father (MoBa). The controls include the mothers of naturally conceived singletons. All the results are shown with the mothers of naturally conceived DZ twins as the reference group. MAR = medically assisted reproductionGroupNP* valueBeta (95% CI)Mean maternal ageMean heightMean BMI before pregnancyNumber of older children% Smoking yesDZ naturally658--31.34168.924.661.9077.46Controls634221.20 × 10^−4^ − 6.33 × 10^7^ (− 9.61 × 10^7^ to − 3.05 × 10^7^)30.12168.3223.122.0981.27DZ MAR3287.70 × 10^−4^ − 1.10 × 10^8^ (− 1.76 × 10^8^ to − 1.76 × 10^7^)32.42168.824.101.4376.06MZ naturally3130.019 − 7.58 × 10^7^ (− 1.43 × 10^8^ to − 1.40 × 10^6^)29.83167.623.951.7380.27Singletons MAR13886.60 × 10^−4^ − 8.63 × 10^7^ (− 1.37 × 10^8^ to − 3.57 × 10^7^)33.17168.6723.361.5678.10Italic inidcates a significant effect based on the alpha of 0.01

Mothers of naturally conceived DZ and MZ twins

In NTR, we observed a significantly higher mean PGS in the mothers of natural DZ twins than in the mothers of natural MZ twins, for whom their twin pregnancy is not the result of a double ovulation (P = 2.09 × 10^−10^; Table 1 and S2; Fig. 3C). In MoBa, the mean PGS is also higher in the mothers of natural DZ twins compared to the mothers of naturally conceived MZ twins based on a nominal p-value of 0.05, but this was attenuated after correcting for multiple testing using the Bonferroni correction as described in the methods (P = 0.019; Table 2 and S3; Fig. 3D). When comparing controls to mothers of natural MZ twins, we observed no significant difference in the PGS for DZT in the NTR (P = 0.350; Table S4) or MoBa (P = 0.84; Table S5; Figure S3 A). In NTR, mothers of MAR MZ twins (P = 0.185; Table S4; Figure S3B) did not differ from the controls.

Mothers of naturally conceived DZ twins and of DZ MAR twins

Mothers of natural DZ twins have a higher PGS on average compared to mothers of MAR DZ twins (P = 1.33 × 10^−6^; Table 1; Fig. 4A). The significant difference between mothers of natural DZ twins and the ART DZ twin mothers was also observed in MoBa (P = 7.70 × 10^−4^; Table 2; Fig. 4B). We then separated the DZ MAR mothers in NTR into a group of mothers who received OI treatment and the IVF/ICSI mothers. Here, we observed that the mothers who received OI treatment, where couple infertility likely has a female component, had a lower PGS for DZT compared to the mothers of naturally conceived DZ twins (P = 2.00 × 10^−9^; Table 1; Fig. 5). The PGS in the other groups of DZ MAR mothers, where couple infertility can be due to male, female, or combined infertility, did not differ from those in the mothers of naturally conceived DZ twins (Table 1), indicating that the difference in PGS between the mothers of naturally conceived DZ twins and the MAR DZ mothers is driven by maternal infertility issues.Fig. 4. The mean PGS for dizygotic (DZ) twinning in mothers of naturally conceived DZ twins and mothers from DZ twins that received medically assisted reproduction (MAR) treatments from A the Netherlands Twin Register (NTR) and B the Norwegian Mother, Children, and Father Cohort Study (MoBa)Fig. 5. The mean PGS of the DZ medically assisted reproduction (MAR) group of the NTR is further split between ovulation induction (OI) and IVF and ICSI

Singleton mothers

No difference was observed between mothers of the natural singletons and MAR singletons in MoBa (P = 0.32; Table S5; Figure S4). The inclusion of covariates, i.e., maternal age, age at first birth, height, BMI before pregnancy, number of older children, and smoking, did not change these results (Tables S2 and S3).

Time to pregnancy

In the mothers of naturally conceived twins (NTR) and singletons (MoBa), we performed regression analyses between the DZT PGS and time to pregnancy measured in months. We found no effect in the mothers of naturally conceived twins (Table S6); however, we observed that the DZT PGS was associated with a shorten TTP (P = 0.0036; Table S7) when looking at TTP as a continuous variable in the much larger group of singleton mothers. In a categorial, per month, analysis in MoBa, we observed a significant decrease of the PGSs for the mothers for whom it took 12 or more months to conceive compared to mothers who had an unplanned pregnancy (P = 0.0138; Table S7).

Genetic correlations

We calculated genetic correlations between DZT and seven fertility related phenotypes including PCOS, anovulatory infertility, endometriosis, number of children ever born (NEB), childlessness (CL), and male infertility (Table 3). DZT correlated negatively with anovulatory infertility (rg = − 0.698; P = 4.93 × 10^−5^) and PCOS (rg = − 0.278; P = 0.0353) and positively with endometriosis (rg = 0.279; P = 0.0033). Table 3. Genetic correlations between dizygotic twinning (DZT) and seven fertility related phenotypes. The reference indicates the original publication of the GWAS. rg = genetic correlationTraitNSNP-based heritabilityP heritabilityr_gSE rgP rg_ReferenceDZT716.517 (8.265 MODZT, 26.252 DZT)0.50%0.001---Mbarek et al. (2024) [6]Anovulatory infertility323.473 (7.359 cases)0.80%0.002 − 0.6980.1625.099 × 10^−5^Venkatesh et al. (2024) [40]Polycystic ovarian syndrome113.248 (10.074 cases)12.40%0.020 − 0.2780.1320.035Day et al. (2018) [39]Miscarriages428.523 (69.054 cases)0.43%0.001 − 0.1650.1510.272Laisk et al. (2020) [42]Endometriosis474.160 (60.674 cases)1.63%0.0020.2790.0950.003Rahmioglu et al. (2023) [41]Number of children ever born785.6042.11%0.0010.0930.0640.147Mathieson et al. (2023) [38]Childlessness450.0823.62%0.0020.0140.0730.853Mathieson et al. (2023) [38]Male infertility3.225 (1.274 cases)11.5%0.3270.0480.1210.694Cerván-Martín et al. (2022) [43]

Discussion

Our results indicate DZT to be a valuable phenotype for studying female fertility and infertility as the DZT PGS distinguishes between twin mothers of naturally conceived pregnancies, and twin mothers who received MAR treatments and can identify those at an increased risk of requiring MAR treatments. The PGS is also associated with a shorter time to pregnancy and is lowest in mothers who received hormonal induction of ovulation, suggesting that the genetic variants associating with DZT primarily affect ovulation processes. This suggestion is supported by the strong genetic correlation between DZT and anovulatory infertility.

We observe that, in contrast to many other human traits, the best performing DZT PGS does not include all possible SNPs [44]. Here, we include the f = 0.001 of the SNPs with a nonzero effect, indicating that DZT is not as poly genic or omnigenic as other complex traits [45]. The GENESIS analysis further confirms that DZT is more likely associated with a smaller number of variants with larger effect sizes as GENESIS predicts 2244 susceptible SNPs for DZT, while it has predicted over 30,000 susceptible SNPS for highly polygenic traits such as major depressive disorder [23]. SbayesS predicts that approximately 0.2% of all SNPs will have a nonzero effect on DZT, while for highly polygenic traits such as educational attainment, more than 5% of common SNPs are predicted to have a nonzero effects [24]. The relatively low polygenicity for DZT found here is in line with findings for other fertility traits in a recent large infertility GWAS [40, 46]. These results imply that the few variants associated with DZT are likely to have a large effect on overall female fertility.

The PGS based on the first, 2016 DZT GWAS, which obtained evidence for FSHB and SMAD3, was associated with having children, greater lifetime parity, and earlier age at first child in the Icelandic population [8]. Here, we add that the PGS for DZT, based on the 2024 DZT GWAS, is higher in mothers who naturally conceived DZ twins compared to DZ twin mothers who received MAR treatments. In the mothers of singletons, the PGS is associated with a reduced time to pregnancy. These findings highlight the connection between DZT and female fertility, emphasizing their potential significance for women contemplating delaying childbirth. While a lower PGS could in future serve as an important indicator for those considering egg freezing at a younger age, it is important to note that this approach may not yet fully feasible for implementation for this purpose.

The genetic signal detected for DZT is most likely driven by genetic influences on ovulation rate, as the largest difference in PGS was observed between the mothers of naturally conceived DZ twins and the mothers who received OI treatment only, which generally, but not exclusively, indicates maternal fertility issues due to anovulation problems. OI is typically the first-line treatment for infertility, and when successful, it therefore suggests ovulatory dysfunction. However, if further treatments such as IVF or ICSI are needed, it usually indicates the presence of additional fertility challenges beyond anovulation [18, 47]. Anovulatory infertility is in direct contrast to the hyper-ovulation seen in mothers of naturally conceived DZ twins, which is also confirmed with the strong negative genetic correlation between DZT and anovulatory infertility estimated at rg = − 0.698 [48, 49]. These findings point to ovulation rate as a genetic spectrum with anovulatory infertility on one end and hyper-ovulation as indicated by natural DZT on the other end.

Pregnancies after ART, especially twin pregnancies, carry higher mortality and morbidity risks for mother and child and more miscarriages are seen with MAR treatments [50–52]. Also, MAR treatments have on average a low success rate of approximately 20% varying between 5 and 40% depending on, e.g., maternal age [19, 53]. Our results support the idea that genetic investigation with respect to DZT could aid in improving the outcomes of ART treatments [54]. For example, combining the DZT PGS with clinical features may serve in better predictions of who will need MAR treatments and which specific treatments has the best chance of success [55].

Female infertility is an overarching term that covers heterogenous problems. DZT was negatively correlated with anovulatory infertility and PCOS (rg = − 0.698 and rg = − 0.278 respectively), but positively with endometriosis (rg = 0.279). The positive correlation with endometriosis seems puzzling, but is in line with a recent study that compared pregnancy outcomes among women with endometriosis and uterine cysts and controls, which observed that mothers with endometrioses have an increased risk for twinning, regardless of MAR status [56].

It has been proposed that female fertility traits are mediated by a common genetic factor, though the direction of effects of these factors seem to vary per fertility trait [57]. A potential common genetic factor is the top locus of the DZT GWAMA, FSHB as this locus has also been associated with endometriosis, PCOS, ovarian cyst, and menorrhagia, and the direction of effect differs between traits [6, 39, 41, 57–59].

Some limitations should be considered. The underlying reasons for ART MAR treatments in NTR and MoBa could be caused by maternal or paternal fertility problems, or a combination of both. Treatment type, which was known in the NTR but not in MoBa, can serve as a proxy for the origin of the fertility problems. For example, as OI most likely indicates maternal problems, but we cannot exclude the possibility that reduced male fertility also plays a role. The lack of further differentiations may also explain why no difference is observed in singleton mothers of the MoBa cohort. Still, in singleton mothers, an increased DZT PGS was associated with a shorter TTP, further supporting the relation between DZT and heightened female fertility.

A strength of our study is that we analyzed data from cohorts from two different countries, including over 7.000 participants in NTR and over 65.000 participants in MoBa with multiple control groups. Still, our sample is solely based on individuals of European ancestries. DZ twinning is a trait with huge global differences in prevalence [60, 61]. To test if the effect is also present in other ancestries, we calculated the PGS in three groups of twin mothers in non-European ethnic groups (mostly Caribbean ancestry; N = 118) from the NTR. Similar to the main findings based on the sample from European ancestry, Figure S5 suggests that the PGS is lower in the DZ twin mothers who received MAR treatments compared to the mothers of naturally conceived DZ twins. However, follow-up based on larger sample sizes is clearly required. GWAS and PGS analyses in diverse samples with high DZ twinning rates, e.g., as observed in Nigeria, can enable tests of genetic explanations for the large global differences.

In summary, we find evidence for a genetic spectrum ranging from anovulation/infertility to becoming the mother of spontaneous DZT. A person’s liability on this scale informs on DZT risk and other female fertility traits such as time to pregnancy, PCOS, and endometriosis and implies that DZT is a valuable model for understanding female infertility.

Supplementary Information

Below is the link to the electronic supplementary material.Supplementary file1 (DOCX 1333 KB)