Genome-wide association study meta-analysis of dizygotic twinning illuminates genetic regulation of female fecundity

Abstract

Study question: Which genetic factors regulate female propensity for giving birth to spontaneous dizygotic (DZ) twins?

Summary answer: We identified four new loci, GNRH1, FSHR, ZFPM1, and IPO8, in addition to previously identified loci, FSHB and SMAD3.

What is known already: The propensity to give birth to DZ twins runs in families. Earlier, we reported that FSHB and SMAD3 as associated with DZ twinning and female fertility measures.

Study design, size, duration: We conducted a genome-wide association meta-analysis (GWAMA) of mothers of spontaneous dizygotic (DZ) twins (8265 cases, 264 567 controls) and of independent DZ twin offspring (26 252 cases, 417 433 controls).

Participants/materials, setting, methods: Over 700 000 mothers of DZ twins, twin individuals and singletons from large cohorts in Australia/New Zealand, Europe, and the USA were carefully screened to exclude twins born after use of ARTs. Genetic association analyses by cohort were followed by meta-analysis, phenome wide association studies (PheWAS), in silico and in vivo annotations, and Zebrafish functional validation.

Main results and the role of chance: This study enlarges the sample size considerably from previous efforts, finding four genome-wide significant loci, including two novel signals and a further two novel genes that are implicated by gene level enrichment analyses. The novel loci, GNRH1 and FSHR, have well-established roles in female reproduction whereas ZFPM1 and IPO8 have not previously been implicated in female fertility. We found significant genetic correlations with multiple aspects of female reproduction and body size as well as evidence for significant selection against DZ twinning during human evolution. The 26 top single nucleotide polymorphisms (SNPs) from our GWAMA in European-origin participants weakly predicted the crude twinning rates in 47 non-European populations (r = 0.23 between risk score and population prevalence, s.e. 0.11, 1-tail P = 0.058) indicating that genome-wide association studies (GWAS) are needed in African and Asian populations to explore the causes of their respectively high and low DZ twinning rates. In vivo functional tests in zebrafish for IPO8 validated its essential role in female, but not male, fertility. In most regions, risk SNPs linked to known expression quantitative trait loci (eQTLs). Top SNPs were associated with in vivo reproductive hormone levels with the top pathways including hormone ligand binding receptors and the ovulation cycle.

Large scale data: The full DZT GWAS summary statistics will made available after publication through the GWAS catalog (https://www.ebi.ac.uk/gwas/).

Limitations, reasons for caution: Our study only included European ancestry cohorts. Inclusion of data from Africa (with the highest twining rate) and Asia (with the lowest rate) would illuminate further the biology of twinning and female fertility.

Wider implications of the findings: About one in 40 babies born in the world is a twin and there is much speculation on why twinning runs in families. We hope our results will inform investigations of ovarian response in new and existing ARTs and the causes of female infertility.

Study funding/competing interest(s): Support for the Netherlands Twin Register came from the Netherlands Organization for Scientific Research (NWO) and The Netherlands Organization for Health Research and Development (ZonMW) grants, 904-61-193, 480-04-004, 400-05-717, Addiction-31160008, 911-09-032, Biobanking and Biomolecular Resources Research Infrastructure (BBMRI.NL, 184.021.007), Royal Netherlands Academy of Science Professor Award (PAH/6635) to DIB, European Research Council (ERC-230374), Rutgers University Cell and DNA Repository (NIMH U24 MH068457-06), the Avera Institute, Sioux Falls, South Dakota (USA) and the National Institutes of Health (NIH R01 HD042157-01A1) and the Genetic Association Information Network (GAIN) of the Foundation for the National Institutes of Health and Grand Opportunity grants 1RC2 MH089951. The QIMR Berghofer Medical Research Institute (QIMR) study was supported by grants from the National Health and Medical Research Council (NHMRC) of Australia (241944, 339462, 389927, 389875, 389891, 389892, 389938, 443036, 442915, 442981, 496610, 496739, 552485, 552498, 1050208, 1075175). L.Y. is funded by Australian Research Council (Grant number DE200100425). The Minnesota Center for Twin and Family Research (MCTFR) was supported in part by USPHS Grants from the National Institute on Alcohol Abuse and Alcoholism (AA09367 and AA11886) and the National Institute on Drug Abuse (DA05147, DA13240, and DA024417). The Women's Genome Health Study (WGHS) was funded by the National Heart, Lung, and Blood Institute (HL043851 and HL080467) and the National Cancer Institute (CA047988 and UM1CA182913), with support for genotyping provided by Amgen. Data collection in the Finnish Twin Registry has been supported by the Wellcome Trust Sanger Institute, the Broad Institute, ENGAGE-European Network for Genetic and Genomic Epidemiology, FP7-HEALTH-F4-2007, grant agreement number 201413, National Institute of Alcohol Abuse and Alcoholism (grants AA-12502, AA-00145, AA-09203, AA15416, and K02AA018755) and the Academy of Finland (grants 100499, 205585, 118555, 141054, 264146, 308248, 312073 and 336823 to J. Kaprio). TwinsUK is funded by the Wellcome Trust, Medical Research Council, Versus Arthritis, European Union Horizon 2020, Chronic Disease Research Foundation (CDRF), Zoe Ltd and the National Institute for Health Research (NIHR) Clinical Research Network (CRN) and Biomedical Research Centre based at Guy's and St Thomas' NHS Foundation Trust in partnership with King's College London. For NESDA, funding was obtained from the Netherlands Organization for Scientific Research (Geestkracht program grant 10000-1002), the Center for Medical Systems Biology (CSMB, NVVO Genomics), Biobanking and Biomolecular Resources Research Infrastructure (BBMRI-NL), VU University's Institutes for Health and Care Research (EMGO+) and Neuroscience Campus Amsterdam, University Medical Center Groningen, Leiden University Medical Center, National Institutes of Health (NIH, ROI D0042157-01A, MH081802, Grand Opportunity grants 1 RC2 Ml-1089951 and IRC2 MH089995). Part of the genotyping and analyses were funded by the Genetic Association Information Network (GAIN) of the Foundation for the National Institutes of Health. Computing was supported by BiG Grid, the Dutch e-Science Grid, which is financially supported by NWO. Work in the Del Bene lab was supported by the Programme Investissements d'Avenir IHU FOReSIGHT (ANR-18-IAHU-01). C.R. was supported by an EU Horizon 2020 Marie Skłodowska-Curie Action fellowship (H2020-MSCA-IF-2014 #661527). H.S. and K.S. are employees of deCODE Genetics/Amgen. The other authors declare no competing financial interests.

Trial registration number: N/A.

Keywords: dizygotic twinning; female fecundity; fertility; genetics; genome-wide association analysis.

Figure 1.

Association results of the genome wide association study (GWAS) meta-analysis of dizygotic (DZ) twinning versus controls. (A) Manhattan plot of single nucleotide polymorphisms (SNP) associations. The -log₁₀ p value (y-axis) for the SNPs are plotted against their physical chromosomal position (x-axis). The dashed red line represents the genome-wide level of significance (P value< 5 × 10⁻⁸). The rs numbers and gene name indicate the chromosomal region attaining genome-wide significance. For plotting purposes, overlapping data points are not drawn for filtered SNPs with a P-value ≥1 × 10⁻⁵. (B) Manhattan plot of gene-based association results as computed by MAGMA. The -log₁₀P value (y-axis) for the genes are plotted against their physical chromosomal position (x-axis). Input SNPs were mapped to 19171 protein coding genes. Genome-wide significance (red dashed line in the plot) was defined at P = 0.05/19171 = 2.61 × 10⁻⁶.

Figure 2.

Regional association plots for genome wide association study (GWAS) meta-analysis, for newly-reported associations. Regional association plots for the four new regions reported in Table 2. The reference single nucleotide polymorphism (SNP) is the most highly associated one reported. The suggested gene in the text is highlighted in red. Loci are: (a) 16q24.2; (b) 8p21.2; (c) 12p11.21; (d) 2p16.3. Approximate conditional analysis in GCTA on lead SNPs (not shown) shows only one significant association in each case. Plots were prepared using a stand-alone copy of LocusZoom 1.3 (http://locuszoom.org/) for the 1000G Release 3 European LD reference, and hg19/Build 37.

Figure 3.

Genetic relationships between the genome wide association study (GWAS) meta-analysis results of dizygotic (DZ) twinning and 28 other traits. Single nucleotide polymorphism (SNP)-based genetic correlations (r_g) were estimated with linkage disequilibrium (LD) score regression. The bars represent the standard error of r_g. The genetic correlation estimates are colour-coded by trait category. BMI, body mass index; SHBG, sex hormone binding globulin; HDL, high-density lipoprotein; ECG, electrocardiogram; ADHD, attention deficit hyperactivity disorder.

Figure 4.

Association of top dizygotic (DZ) twinning single nucleotide polymorphisms (SNPs) with in vivo reproductive hormone data from four independent studies. Heat map showing the beta value of the association of each DZ twinning SNP and corresponding hormone data supplied by the indicated study. SNPs with a positive effect (positive beta) are shown in red, SNPs with a negative effect (negative beta) are shown in blue. Darker shading indicates beta values of larger magnitude. A single asterisk (*) indicates nominal significance at alpha = 0.05, and the double asterisk (**) indicates significance at the alpha corrected for multiple testing, i.e. alpha = 0.05/(10 × 17). Cells in light gray represent missing/unavailable data. Additional details, including study sample sizes and hormone measurement specifications, can be found in Supplementary Table S12. LHCGH, luteinizing hormone/choriogonadotropin hormone; FSHR, follicle stimulating hormone receptor; HLA-G, human leukocyte antigen G; GNRH1, gonadotropin releasing hormone 1; FSHB, follicle stimulating hormone subunit beta; IPO8, importin 8; SMAD3, mothers against decapentaplegic homolog 3; ZFPM1, zinc finger protein; FOG, Family Member 1; SHBG, sex hormone binding globulin; FSH, follicle stimulating hormone; LH, luteinizing hormone; Test, total testosterone; Est, estradiol. Iceland (deCODE), NESDA, TwinsUK and UKB are the four studies.

Figure 5.

The correlation between population twinning rates from 53 populations and predicted dizygotic twin (DZT) polygenic risk score (PRS) for each population based on European effect sizes of top single nucleotide polymorphisms (SNPs) from Table 2. The twinning rate data are derived from Smits & Monden Demographic and Health Surveys (Smits and Monden, 2011). The HGDP group twinning rates are taken as those of the country from which they derive (e.g. Mandinka from Senegal); we have added twinning rate data for UK and Japan from Imaizumi (Imaizumi, 2003), other European twinning rates from Heino et al. (2016) and rates for Australia from the Australian Bureau of Statistics, 2010. The SNPs and their effect sizes are the 26 (unique) top associated SNPs from our genome wide association meta-analysis (GWAMA) listed in Table 2, excluding (for technical reasons) rs17293443, rs16050687, and rs17293443. Allele frequencies were obtained from the Human Genome Diversity Project and the 1000 Genomes Projects. We calculated the Pearson correlation and a bootstrap standard error (R boot package, 1000 replicates) using country-level rate and mean PRS for 53 populations. The vertical dotted line is the mean PRS for all mothers of DZ twins in the current study. The oblique line shows the regression line of observed twinning rate on the predicted DZT PRS for each population.

Figure 6.

ipo8 deficiency in zebrafish compromises maternal, but not paternal, fertility. (A) Representative wild type (left) and maternal zygotic ipo8 knockout (right) zebrafish embryo at 2.5 hpf. Embryos devoid of ipo8 show a less compact dome with loosened blastomeres, as indicated by the black arrows. Scale bar 250 µm. (B) At the 64-cell stage, prior the onset of zygotic transcription, endogenous ipo8 mRNA, revealed by whole-mount in situ hybridization, is contributed through maternal deposition in the egg and is absent in embryos derived from homozygous ipo8 knockout females. Scale bar 250 µm. (C) Endogenous ipo8 mRNA present in the 64-cell stage, measured by QPCR, originates from the maternal, but not paternal, contribution. (D) The size of clutches laid by wildtype, heterozygous or homozygous ipo8 knockout mothers is not statistically different. Paternal ipo8 genotype has no effect on egg clutch size. N > 6 for each cross. (E) Survival of embryos at 24 h post fertilization is solely dependent on maternal, but not paternal, presence of wild type ipo8. All error bars show 95% confidence interval. ∗∗P = <0.01, ∗∗∗P = <0.001; Mann–Whitney test.