Spontaneous dizygotic (DZ) twinning rates vary markedly between the major ethnicities (Bulmer, Reference Bulmer1970; Monden et al., Reference Monden, Pison and Smits2021; Sear et al., Reference Sear, Lawson, Kaplan and Shenk2016). The DZ twinning rate can be a useful index of fertility in a population (Tong & Short, Reference Tong and Short1998), and age-specific changes in DZ twinning may reflect declining ovarian follicle reserve and individual fertility (Beemsterboer et al., Reference Beemsterboer, Homburg, Gorter, Schats, Hompes and Lambalk2006). On the other hand, twinning is associated with increased risks to maternal and infant health (Monden & Smits, Reference Monden and Smits2017; Santana et al., Reference Santana, Cecatti, Surita, Silveira, Costa, Souza, Mazhar, Jayaratne, Qureshi, Sousa and Vogel2016; Santana et al., Reference Santana, Silveira, Costa, Souza, Surita, Souza, Mazhar, Jayaratne, Qureshi, Sousa, Vogel and Cecatti2018) and it is important to understand the factors regulating DZ twinning frequency and their relationship to fertility and reproductive aging. We previously reported the first two significant loci for being a mother of spontaneous DZ twins (MoDZT), based on a GWAMA of 1980 cases plus 12,953 controls (Mbarek et al., Reference Mbarek, Steinberg, Nyholt, Gordon, Miller, McRae, Hottenga, Day, Willemsen, de Geus, Davies, Martin, Penninx, Jansen, McAloney, Vink, Kaprio, Plomin, Spector and Boomsma2016). We found loci close to two genes, FSHB, the structural locus for follicle stimulating hormone (FSH) beta subunit and SMAD3, which regulates the response of the ovaries to FSH, and replicated these findings in gene-based testing of ‘being a DZ twin’ in the UK Biobank. We have now expanded the Australian and New Zealand component of that study to 3273 MoDZT plus 24,009 controls of European ancestry.
Methods
The analysis was a case/control genome-wide association study (GWAS), where the cases were MoDZT and ‘controls’ were those not known to be cases and who have no self-reported closely related DZ twins (based on known pedigrees, or relatedness with cases above pi-hat ∼0.1 in an Identity By Descent (IBD) analysis; see later description under ‘Controls’). No phenotype information was used apart from zygosity of the twins, assisted reproductive technology (ART) status (see below) and family history of DZ twinning, as foregoing. All studies were approved by the QIMR Berghofer Human Research Ethics Committee.
Cases
Cases were genotyped MoDZT drawn from a twin family dataset collected and held by QIMR Berghofer (‘QIMR’). This primarily included (a) families with twins from various studies that had recruited twins and members of their families since 1978 via the Australian Twin Registry and from whom DNA samples were subsequently collected (most of these twins were born before 1972 so ART is not an issue); (b) families with twins born after 1980 and enrolled in the Queensland Twin Registry (the Brisbane Adolescent Twin sample (BATS), Project 5 in Medland, Nyholt et al. (Reference Medland, Nyholt, Painter, McEvoy, McRae, Zhu, Gordon, Ferreira, Wright, Henders, Campbell, Duffy, Hansell, Macgregor, Slutske, Heath, Montgomery and Martin2009); see also Wright and Martin (Reference Wright and Martin2004) and Medland, Duffy et al. (Reference Medland, Duffy, Wright, Geffen, Hay, Levy, van-Beijsterveldt, Willemsen, Townsend, White, Hewitt, Mackey, Bailey, Slutske, Nyholt, Treloar, Martin and Boomsma2009); (c) the QIMR Twinning Genetics (TG) study, which recruited MoDZT from multiplex families with closely related MoDZT, mainly sisters whom we had previously used in a sib-pair linkage study of DZ twinning. Participants were from Australia or New Zealand and ascertained through Mothers-of-Twins clubs and elsewhere (Painter et al., Reference Painter, Willemsen, Nyholt, Hoekstra, Duffy, Henders, Wallace, Healey, Cannon-Albright, Skolnick, Martin, Boomsma and Montgomery2010).
The two big issues with acceptability of a MoDZT case for our study are (1) that her twins are definitely DZ, and (2) that she conceived the twins spontaneously and did not use any ART. Both issues are addressed below. A schema summarizing our ascertainment procedure with regard to both ART and zygosity is shown in Appendix A.
Controls
Controls were drawn from two groups from the general Australian population, also filtered genetically to include only Europeans: (1) individuals from the studies supplying case groups (a) and (b) above where there were no known DZ twins in the family and genotyping was available; (2) genotyped individuals from the QIMR Berghofer QSkin Study (Olsen et al., Reference Olsen, Green, Neale, Webb, Cicero, Jackman, O’Brien, Perry, Ranieri and Whiteman2012) which is an unselected sample of Queenslanders aged 40−70 randomly contacted via the electoral roll (voter registration is compulsory in Australia and compliance is 97%). Around 44,000 completed an online survey on various aspects of health, but with a focus on skin cancer. All were invited to donate a saliva sample. QSkin subjects who were genotyped and had consented to re-contact were contacted by email and asked to complete a brief online screening questionnaire (Appendix B) about their family history of twinning including (female participants) whether they themselves had had twins and if so, their ART status and zygosity of the twins. Of the genotyped 17,642 QSkin subjects, we excluded 378 due to non-European ancestry, followed by 897 who were DZ twins themselves, a mother of twins, or had a close twin relative (leaving 16,367 at this point). For QSkin and elsewhere, we retained controls even if they were male, after running the final analysis twice: first with both male and female controls (N = 24,009, the main analysis), and second with only female controls (N = 12,819 controls). This showed that the resulting near doubling of sample size for controls outweighs the fact that male controls, by virtue of their inability to become pregnant, were only screened for a predisposition towards DZ twinning, based on limited family history information (see Results).
In the final GWAS analysis, tight relatedness exclusions were applied to controls who had an IBD-based pi-hat ≥ 0.1 with any DZ twin or case/mother of DZ twin (even unselected cases) in a GRM including all individuals with available genotyping. In these instances the controls were removed, as they greatly outnumbered cases. Genotyped controls (both sexes) remaining after all QC, comprised QSkin (N = 14,577) and other controls as in (1) above, of which N = 658 were from the same batch of GSA genotyping as QSkin, N = 3209 typed on Core+Exome, PsychArray, Omni2.5, OmniExpress chips; N = 2442 typed on the 610K chip, and N = 3123 typed on 317K or 370K chips for a total of total N = 24,009 controls (12,820 female, 11,189 male).
Zygosity – Ensuring Twins Were DZ
For opposite-sex pairs (about one third of spontaneous European twins), zygosity is unequivocally DZ, but for same-sex pairs, zygosity was initially based on responses to standard questions regarding similarity of appearance. Later this was refined by genetic testing, initially using the Profiler© microsatellite set but more recently using genomewide genotyping with a SNP array so that in most cases final zygosity is based on direct genetic evidence. Direct genotyping of the twins and a variety of clinical or laboratory tests performed at and after recruitment already made the zygosity of twins for whom DNA was available quite certain. However, in many cases, the twins of a mother selected as a potential case had not been genotyped so we could not genetically confirm zygosity. In these cases, zygosity was established by questionnaires containing standard zygosity items and, in cases of remaining doubt, by telephone interview. If, after these steps, zygosity was still in doubt then the mother was excluded from the MoDZT case sample.
Screening Out MoDZT Who Had Used ART
There has been a large increase in DZ twinning following the widespread adoption from the 1980s on of ART, such as in vitro fertilization (IVF); before 1980 use of ART was uncommon (Monden et al., Reference Monden, Pison and Smits2021). In our previous work (Mbarek et al., Reference Mbarek, Steinberg, Nyholt, Gordon, Miller, McRae, Hottenga, Day, Willemsen, de Geus, Davies, Martin, Penninx, Jansen, McAloney, Vink, Kaprio, Plomin, Spector and Boomsma2016) we showed that although ART accounts for only a small fraction of twin pregnancies, the sensitivity of the phenotype definition to contamination with DZ twins arising after ART is high, with as few as 10% ART cases entirely ablating any association signals. In addition to dilution of a genetic signal due to phenocopies, it is possible that some women requiring ART to conceive are genetically predisposed to low fertility and their inclusion would therefore tend to cancel the very genetic signals we are looking for.
For mothers contemplated for inclusion as MoDZT cases, where twins were born after 1980 and we had no existing ART screening information, we therefore made a substantial effort to re-contact the mother by telephone or email to confirm their ART status (see questionnaire in Appendix C) and at the same time to confirm twins’ zygosity. Phase 1 of this effort involved attempting to contact 755 mothers of twins in the youngest (mostly QTwin-derived) cohort with 397 (53%) completing screening questions. Phase 2 involved contacting other mothers targeted for new GSA genotyping outside of the original Twinning Genetics Study; with 484 mothers contacted successfully; 81 of whom failed the ART exclusion and/or zygosity check; with another 427 uncontactable. Those failing or uncontactable were excluded from further genotyping, with an initial highest-priority candidate list of 793 individuals being reduced by 78 (9.8%) due to screening questions and 177 (22.3%) by lack of contact (538 remaining). An overall list of known ART cases from both phases was also excluded from the analysis itself. The same questionnaire also confirmed that pre-1980 use of ART was indeed rare, which gives some confidence in assuming lack of ART in that age group.
These efforts, along with existing screening data for other individuals, will have removed the great majority of potential cases where twin pregnancies were due to ART or twins were not DZ. However, there remain some mothers with pre-existing genotyping who could not be contacted and have not been screened for ART or re-contacted for zygosity. We chose to retain these uncontactable cases and so not filter them out explicitly, on the basis that (1) the great majority of these births were pre-1980 and unlikely to use ART; (2) zygosity error is known to be quite low.
Genotyping
These data are held as an integrated dataset (R10) expanded (added 660K, Omni family and GSA genotyping including QSkin) and reimputed from that used in Medland, Nyholt, et al. (Reference Medland, Nyholt, Painter, McEvoy, McRae, Zhu, Gordon, Ferreira, Wright, Henders, Campbell, Duffy, Hansell, Macgregor, Slutske, Heath, Montgomery and Martin2009); also used in Wood et al. (Reference Wood, Esko, Yang, Vedantam, Pers, Gustafsson, Chu, Estrada, Luan, Kutalik, Amin, Buchkovich, Croteau-Chonka, Day, Duan, Fall, Fehrmann, Ferreira and Jackson2014) and elsewhere. Genotyping was taken as a subset of the imputed version of an existing dataset of ∼51,834 individuals. These were genotyped in batches of up to ∼22,000 on successive generations of Illumina SNP arrays from 2006 onward. Genotyping as a result was spread across three families/generations of Illumina chip (listed chronologically): (1) HapMap-derived 317K/370K/610K; and (2) 1000 Genomes-derived, Omni family (Core+Exome, PsychChip, Omni 2.5, Omni Express). (3) Global Screening Array (GSA v1); this latter includes Qskin control samples which were typed on the Illumina Global Screening Array (model GSAMD-24v1-0_20011747_A1) array at Erasmus University Medical Centre’s Human Genomics Facility (HuGe-F, André Uitterlinden director). It also includes a batch of 1922 samples (1806 after QC), including 775 MoDZT sister pairs, 329 MoDZT from parent-child trios, 309 singleton MoDZTs and 509 MoDZT from other studies, which were genotyped at the Avera Institute of Human Genetics, Sioux Falls SD, at Avera Health using the Avera-NTR GSA array (model GSA-24v1-0_A1_avera_20170513f, a semi-custom SNP array made by Illumina (Beck et al., Reference Beck, Hottenga, Mbarek, Finnicum, Ehli, Hur, Martin, de Geus, Boomsma and Davies2019).
The entire dataset contains N = 19,428; 14,210; 23,821 individuals from HapMap-derived Illumina arrays; 1000G-derived Illumina arrays; and GSA v1 respectively. Analyzed individuals are a subset of this with N = 6322; 3922; 17,038 individuals respectively.
Individual batches of genotyping were called and QC’d as per standard protocols (within and post GenomeStudio, including dropping individuals <97%, and SNPs with GenTrain Score < 0.6, MAF<1%, p(HWE)<10-6, mean(GC) <0.7, call rate <95%) before being integrated into an overall dataset. Further Mendelian and relatedness checks were performed by Identity By Descent (IBD) calculation run on the overall dataset. Batch effects at the genotyping stage are not significant after QC; they were checked by various means including allele frequency comparisons and repeat genotyping of the same samples. Family-based data cleaning prior to GSA genotyping allowed most sample issues to be identified and resolved. Other misidentified or wrong sex samples were removed. The X chromosome data were available post-QC.
Data Cleaning and Imputation
A Principal Components Analysis (PCA) was run using smartpca (in EigenSoft 7.2.0; original version described in Price et al., Reference Price, Patterson, Plenge, Weinblatt, Shadick and Reich2006) with samples from HapMap Phase 3 and Genome-EUTWIN populations used to define the PC axes; all QIMR genotyping was projected onto those axes. About 3% of individuals (from the overall dataset) with PC1 and/or PC2 >6sd from the mean of European reference populations were excluded as ancestry outliers. In the association analysis, PC1-PC4 were also used as covariates in all analyses, although from past experience, there is minimal population stratification after the stated exclusions.
Imputation was run on the Michigan Imputation Server using SHAPEIT (phasing), minimac3 (imputation) and the HRCr1.1 reference panel, for each of the three SNP-array families. Each run used only individuals genotyped with those arrays, and observed markers passing QC in all relevant batches (19,428 individuals, 280,280 markers HapMap-derived; 14,210 individuals, 238,591 markers 1000G-derived; 23,821 individuals, 481,926 markers GSA) and merged post-imputation, preferencing HapMap-derived; 1000G-derived; GSA in that order for each individual. Two binary analysis covariates record that choice of imputation run, to model any differences in imputation between the chip families.
Association Analysis (GWAS)
Although controls closely related to cases (pi-hat ≥0.1) were removed to improve the quality of controls, cases are drawn primarily from family-based studies, as are the non-QSkin controls, so substantial relatedness still exists within both the case and control groups. Reducing the analysis to an ‘unrelateds only’ analysis would have considerably reduced the sample size and power of the analysis. Hence the association tests were run using the R package SAIGE (v 0.36.3.3) (Zhou et al., Reference Zhou, Nielsen, Fritsche, Dey, Gabrielsen, Wolford, LeFaive, VandeHaar, Gagliano, Gifford, Bastarache, Wei, Denny, Lin, Hveem, Kang, Abecasis, Willer and Lee2018), which uses a Generalised Linear Mixed Model (GLMM) to model relatedness as a part of the initial stage of its analysis, and also has a high tolerance to unbalanced case control ratios as applies here. SAIGE ‘Step 1’ (the null GLMM fit) was run with, as input, the phenotype and covariates plus a genomewide set of observed genotypes for ∼39,500 markers available in all batches of genotyping. Phenotype was the case/control variable coded 1 = case (N = 3273), 0 = control (N = 24,009). Covariates were the first 4 genetic Principal Components (PCs) from the smartpca ancestry PCA run (PC1-PC4), plus two binary covariates encoding which of the three imputation runs was used for the individual in question. SAIGE ‘Step 2’ was then run in parallel (for speed) across blocks of ∼50,000 markers with as input the Variance Ratio and GMMAT model estimate files from ‘Step 1’ plus imputed genotypes (against the HRC r1.1 reference panel; or TOPMed r2, females-only, for the X chromosome) for that block of markers. Step 2 used a minor allele frequency (MAF) cutoff of 0.001 and minimum minor allele count (MAC) cutoff of 5. The results for these blocks were then concatenated and filtered based on imputation metadata. For the purposes of results presented here, only markers with MAF ≥ 0.01 (1%) and r 2 ≥ 0.6 (the lowest, that is, worst r 2 across the three imputation runs) were retained.
Results
The results of the SAIGE SNP-based GWAS analysis are summarized in Figure 1 (a Manhattan Plot, i.e., p value vs. genomic position) and Figure 2 (the corresponding Quantile-Quantile or Q-Q Plot), after the chosen filter of MAF ≥ 0.01 (1%) and Rsq ≥ 0.6. These show four unique association peaks below the conventional genome-wide significance threshold for GWAS results, p = 5 × 10-8, along with a number of narrowly below significant associated regions that are not explored further here. The significant regions are listed in Table 1. The Genomic Inflation Factor (λ) is ∼1.125 (1.131 for autosomes). No significant peaks are seen on the X-chromosome.
Figure 1. Manhattan plot for additive model (MoDZT vs. controls). Shows association p values calculated by SAIGE (-log10 transformed) after removing markers with MAF < 1% or Rsq < 0.6. Genomewide significance threshold p = 5 × 10-8 is the upper line. Suggestive threshold 10-5 is the lower line. Horizontal axis is genomic position, measured in basepairs, increasing toward the right (hg19). Labels match the genes quoted for the adjacent peak in Table 1.
Note: MoDZT, mothers of dizygotic twins; MAF, minor allele frequency.
Figure 2. Q-Q plot for additive model (MoDZT vs. controls). Shows association p values calculated by SAIGE (-log10 transformed, vertical axis) after removing markers with MAF < 1% or Rsq < 0.6; against p value expected at that rank. The red line and shaded area mark the 95% confidence interval under the null (λ = 1). Observed λ = 1.131.
Note: MoDZT, mothers of dizygotic twins; MAF, minor allele frequency.
Table 1. Genomewide significant genes and loci for additive model (MoDZT vs. control), after removing markers with MAF < 1% or Rsq < 0.6. Positions are for NCBI Human Genome Build 37 (also known as hg19), as chromosome:basepair_position. Only the most associated SNP for each gene/region is shown; each association peak has up to 221 genomewide significant SNPs (for FSHB) or 19 (for SMAD3 case).
Note: MoDZT, mothers of dizygotic twins; SNP, single nucleotide polymorphism.
Figure 3 explores the question of whether to include male controls in the analysis; cases are, by definition, female, so it might be intuited that one should only use female controls. The equivalent GWAS was run for the autosomes using only female controls (N = 12,819; λ for this female-only analysis is ∼1.085) to compare with the analysis using controls of both sexes (N = 24,009); both analyses used N = 3273 cases. The two sets of (QC’d) p values are plotted against each other and show that the analysis using female controls only consistently results in less significant p values for markers below, or close to, the genomewide significance threshold (5 × 10-8). This is particularly notable for the most-associated FSHB and SMAD3 regions; however, there is no pervasive shift from the y = x diagonal for nonassociated high p markers. We take this as our justification for including male controls in our GWAS analyses to maximize power. This is to be expected for a low prevalence (∼1%) trait like DZ twinning where screening out ‘affecteds’ (which we could not do perfectly) makes little difference.
Figure 3. A comparison plot, marker by marker, between p values for that marker in the main analysis with both male and female controls (horizontal axis; 3273 cases, 24,009 controls) and an alternative analysis with only female controls but otherwise equivalent (vertical axis; 3273 cases, 12,819 controls). The diagonal shows y = x. Labels show indicative locations of the three strongest association peaks in the main analysis.
MAGMA v1.0.9a (de Leeuw et al., Reference de Leeuw, Mooij, Heskes and Posthuma2015) was used to perform gene-based association tests using the supplied 1000 European reference genome for LD information. A Manhattan plot showing the resulting p values for the 19,104 genes returning a result is shown in Figure 4. The equivalent tests were also run using VEGAS v2.2 (Mishra et al., Reference Mishra and Macgregor2015) with broadly consistent results for the most associated genes (not presented here). Both programs produce results that are broadly consistent with the strongest SNP associations (Table 1) although some regions are significant at the gene-based level but not at the individual SNP level.
Figure 4. Manhattan plot for MAGMA gene-based test (showing field P_MULTI – composite gene-based p value. The horizontal line is at p = 2.61 × 10-6 = 0.05/N (N = number of tested genes = 19,154). Labels show the most-associated significant or near-significant gene(s) adjoining the label.
The four per-SNP association peaks are shown as regional association plots, with LD information, as Figure 5 (a−d). For comparison purposes, Figure 5(f) also shows a version of panel (a) referenced to the previously published ‘peak’ SNP rs11031006 for the FSHB peak; and panel (g) shows the remaining ‘peak’ SNP rs12064669 from Table 3 of Mbarek et al. (Reference Mbarek, Steinberg, Nyholt, Gordon, Miller, McRae, Hottenga, Day, Willemsen, de Geus, Davies, Martin, Penninx, Jansen, McAloney, Vink, Kaprio, Plomin, Spector and Boomsma2016), which no longer appears associated and is likely to have been a false signal.
Figure 5. Regional association plots and gene annotation for association peaks. Panels (a)-(d) show the four dominant peaks in Table 1, and panel (e) the chr X peak; reported SNP shown as a purple diamond. Panels f-g shows the additional two peaks reported in the previous meta-analysis paper (Mbarek et al., Reference Mbarek, Steinberg, Nyholt, Gordon, Miller, McRae, Hottenga, Day, Willemsen, de Geus, Davies, Martin, Penninx, Jansen, McAloney, Vink, Kaprio, Plomin, Spector and Boomsma2016); (f) rs11031006/11:30226528, which is a lesser-associated chromosome 11 SNP also visible in panel (a), but was the lead chromosome 11 marker reported in Table 3 of the previous paper); and (g) rs12064669/1:230688643 which was reported in that paper, but is not associated in the current analysis. Panels (h)-(i) show the peak surrounding ADRB2 both in the additive model (h) and dominant model (i), referenced to the SNP most-associated in the additive model (rs4705276). All plots prepared using LocusZoom v1.4 with LD r2 estimates derived from 1000 Genomes v3 European genotypes. Grey indicates lack of LD information versus the peak marker (chromosome X and some individual markers).
Note: SNP, single nucleotide polymorphism.
The strongest two associations are those on chromosomes 11 and 15, which reproduce the strongest two associations from our previous paper (Mbarek et al., Reference Mbarek, Steinberg, Nyholt, Gordon, Miller, McRae, Hottenga, Day, Willemsen, de Geus, Davies, Martin, Penninx, Jansen, McAloney, Vink, Kaprio, Plomin, Spector and Boomsma2016) with the expected higher statistical significance. We previously concluded that both relate to genes FSHB (Follicle Stimulating Hormone Beta subunit) and SMAD3 (SMAD Family Member 3) respectively, which are both strong candidates as DZ twinning genes. The MAGMA gene-based tests strongly support both genes along with (at lower significance) the adjacent genes ARL14EP (ADP Ribosylation Factor Like GTPase 14 Effector Protein) and MPPED3 in the case of FSHB (Table 1).
The most convincing of the two other associated regions is at rs4426338, a novel locus close to the gene ZFPM1 (Zinc Finger Protein, FOG Family Member 1) on chromosome 16. The LD block also covers the smaller genes ZNF469 and MIR5189 (see Figure 5c). In fact, the MAGMA gene-based p value for ZNF469 is slightly lower (2.14 × 10-8 vs. 5.26 × 10-8 for ZFPM1; Table 1).
The remaining genomewide significant association at rs6426385 on chromosome 1 has no other marker (after our QC filter) with a p value below 2.85 × 10-4, far higher than the genomewide threshold 5 × 10-8 (Figure 5d). The lack of other associated markers nearby suggests that this is a false-positive result and it is not considered further here. The only nearby gene is noncoding RNA gene LINC01346.
Also prominent and at or near significance in the gene-based results (Table 1; Figure 4) but not reaching genomewide significance for individual SNPs, are ADRB2 on chromosome 5; DOCK5/GNRH1/KCTD9 on chromosome 8; IPO8/CAPRIN2 on chromosome 12; CLYBL on chromosome 13; ADM5/PRMT1/IRF3 on chromosome 19; and RLIM/SLC16A2 on chromosome X. ADRB2 also features in the Dominance Model results (see below). It and CLYBL are close to SNP-wise genomewide significance (Table 1).
Fitting a dominance model. Bulmer (Reference Bulmer1970) interpreted his mother-daughter recurrence data (tetrachoric correlation r = .08) being significantly less than that of sisters (r = .14, p = .007) as implying a strong nonadditive genetic contribution to DZ twinning. We therefore repeated our GWAS using a dominance model, by rerunning SAIGE for the purpose of this exploratory study, with the input dosage values recoded from the imputed genotype probabilities to treat the heterozygote genotype as a double minor allele (i.e., a minor allele dosage of 2 rather than the normal 1 for heterozygotes). The results were effectively identical to the main additive model except for the region in/around the ADRB2 gene (ADRenoceptor Beta 2) on chr 5, with the p value for lead SNP rs4705276 near this gene changing from 1.04 × 10-7 under an additive model to 2.02 × 10-8 under a dominance model (note: run against TOPMedr2 instead of HRCr1.1 used elsewhere). Interestingly, ADRB2 is also the most associated gene in the MAGMA gene-based results, after the chr 11 peak (gene-based p ∼ 2.02 × 10-8). Refer also to Figure 5h-i.
Discussion
Our results, with almost double the number of MoDZT cases (3273 vs. 1980) over the previous Mbarek et al. (Reference Mbarek, Steinberg, Nyholt, Gordon, Miller, McRae, Hottenga, Day, Willemsen, de Geus, Davies, Martin, Penninx, Jansen, McAloney, Vink, Kaprio, Plomin, Spector and Boomsma2016) article, have consolidated and expanded the list of genes involved in spontaneous DZ twinning, at least in Europeans. However, it should be noted that ∼600 cases in the present study also contributed to the meta-analysis in the 2016 paper, so the results of the two studies are not wholly independent.
For the FSHB peak (Figure 5a), we find a different peak SNP (rs74485684), which is only ∼16 kb away from our previously-reported rs11031006, and well within the same LD block, just ∼10kb upstream of FSHB. FSHB is a strong DZ twinning candidate gene as FSH plays a central role in regulating ovarian follicle growth and ovulation (e.g., Trevisan et al., Reference Trevisan, de Oliveira, Christofolini, Barbosa and Bianco2019). Previously we confirmed in an independent Icelandic population that rs11031006 is associated with changed FSH levels (Mbarek et al., Reference Mbarek, Steinberg, Nyholt, Gordon, Miller, McRae, Hottenga, Day, Willemsen, de Geus, Davies, Martin, Penninx, Jansen, McAloney, Vink, Kaprio, Plomin, Spector and Boomsma2016). However, a role for the other gene under the association peak, ARL14EP, cannot be ruled out.
For the SMAD3 peak (Figure 5b), we reproduce the previously found peak SNP rs17293443, which is within that gene’s boundaries. SMAD3 is suspected (based on mouse studies and experiments on human tissue in culture) to be strongly expressed in human ovarian luteinized granulosa cells; to stimulate production of estradiol and progesterone; with possible roles in regulating FSHR (Follicle Stimulating Hormone Receptor) and other genes (e.g., Li et al., Reference Li, Schang, Boehm, Deng, Graff and Bernard2017; Liu et al., Reference Liu, Chen, Xue, Shen, Shi, Dong, Zhang, Liang, Li and Xu2014).
The lead SNP on chromosome 16 near ZFPM1 (Figure 5c) rs4584807 is in high LD with nearby SNPs associated with age at menopause (Mbarek et al., Reference Mbarek, Steinberg, Nyholt, Gordon, Miller, McRae, Hottenga, Day, Willemsen, de Geus, Davies, Martin, Penninx, Jansen, McAloney, Vink, Kaprio, Plomin, Spector and Boomsma2016; Mbarek et al., Reference Mbarek, van de Weijer, van der Zee, Ip, Beck, Abdellaoui, Ehli, Davies, Baselmans, Nivard, Bartels, de Geus and Boomsma2019), sex hormone binding globulin (Ruth et al., Reference Ruth, Day, Tyrrell, Thompson, Wood, Mahajan, Beaumont, Wittemans, Martin, Busch, Erzurumluoglu, Hollis, O’Mara, McCarthy, Langenberg, Easton, Wareham, Burgess, Murray, Ong, Frayling and Perry2020), and FSH concentrations (Chen et al., Reference Chen, Tao, Gao, Zhang, Hu, Mo, Kim, Yang, Tan, Zhang, Qin, Li, Wu, Zhang, Zheng, Xu, Mo and Sun2013; Mbarek et al., Reference Mbarek, Steinberg, Nyholt, Gordon, Miller, McRae, Hottenga, Day, Willemsen, de Geus, Davies, Martin, Penninx, Jansen, McAloney, Vink, Kaprio, Plomin, Spector and Boomsma2016).
The chromosome X gene-based association peak is near RLIM (see also Figure 5e), known also as RNF12, which was first identified as an X-linked activator of X chromosome inactivation (Jonkers et al., Reference Jonkers, Barakat, Achame, Monkhorst, Kenter, Rentmeester, Grosveld, Grootegoed and Gribnau2009). This E3 ubiquitin ligase is involved in the regulation of LIM-homeodomain transcriptions factors. It has been linked to a variety of biological activities, including the appropriate expression of GnRH receptor (GnRHR), which is necessary for the correct regulation of the gonadotropins, LH and FSH, by GnRH (Bach et al., Reference Bach, Rodriguez-Esteban, Carrière, Bhushan, Krones, Rose, Glass, Andersen, Izpisúa Belmonte and Rosenfeld1999; McGillivray et al., Reference McGillivray, Bailey, Ramezani, Kirkwood and Mellon2005). Also nearby is gene SLC16A2, which is involved in transport of thyroid hormones.
Fitting a dominance model, we also find significant evidence for a fourth gene, ADRB2 (Figure 5h-i), which is a plausible twinning candidate. ADRB2 belongs to the family of the beta-adrenergic receptor. It has been shown that both alpha and beta-adrenergic receptors play a role in rodent ovulation (Kannisto et al., Reference Kannisto, Owman and Walles1985, p. 361: ‘The experiments indicate that both a- and ß-adrenergic receptors mediate an increased rate of ovulation through an effect, on the one hand, at the level of the follicle wall and, on the other hand, by a humoral-type of ovarian mechanism.’). Previously, evidence was found pointing to a role for ADRB2 specifically in rat ovulation (Ratner et al., Reference Ratner, Weiss and Sanborn1980). ADRB2 is expressed in human and monkey ovaries and might stimulate both growth of small follicles and induce FSHR, contributing to follicular development (Föhr et al., Reference Föhr, Mayerhofer, Sterzik, Rudolf, Rosenbusch and Gratzl1993; Mayerhofer, et al., Reference Mayerhofer, Dissen, Costa and Ojeda1997; Merz et al., Reference Merz, Saller, Kunz, Xu, Yeoman, Ting, Lawson, Stouffer, Hennebold, Pau, Dissen, Ojeda, Zelinski and Mayerhofe2015).
Our article also addresses the perennial question of whether, for a sex-limited trait like DZT, it is better to use same-sex only controls, or controls of both sexes. In accord with expectations for a low prevalence trait, our analyses show unequivocally that we obtain more power using both sexes as controls, even though males cannot be screened for the trait.
Our results are presented here in bare outline and only briefly, because our immediate intention is to contribute them to a large new meta-analysis, to be published shortly (Mbarek et al., Reference Mbarek, Gordon, Mortlock and Duffyin press).
Acknowledgments
The authors wish to acknowledge the valuable contribution of Joy Brown, a community volunteer (and a mother of DZ twins herself) who carried out most of the work identifying and recruiting mothers of twins from New Zealand, in her own time.
Author contribution
NGM, DLD, GWM designed the study; NGM, SDG, SEM, GWM analysed data; NGM, SDG, HM drafted the manuscript; NGM and JB contributed to recruitment; DCW, CMO, KM, JMA, NAG, SMC, SEL, JJB, GWM contributed to data and sample collection and/or processing.
Financial support
None.
Competing interests
None.
Ethical standard
The authors assert that all procedures contributing to this work comply with the ethical standards of the relevant national and institutional committees on human experimentation and with the Helsinki Declaration of 1975, as revised in 2008.
Appendix A
Summary of criteria for case and control selection
Appendix B
The text of the screening self-report questionnaire for QSkin controls (reformatted)
Appendix C
The text of the screening self-report questionnaire for non-QSkin candidate MoDZTs. Contact was attempted where the birth occurred after 1980 and there was no existing ART information; and/or zygosity was uncertain; for example, ungenotyped same-sex twins.
Open access
