Identification of common genetic variants influencing spontaneous dizygotic twinning and female fertility Hamdi Mbarek,1,21,* Stacy Steinberg,2,22 Dale R. Nyholt,3,4,22 Scott D. Gordon,4 Michael B. Miller,5 Allan F. McRae,6,4 Jouke Jan Hottenga,1,21 Felix R. Day,20 Gonneke Willemsen,1,21 Eco J. de Geus,1,21 Gareth E. Davies,7 Hilary C. Martin,8 Brenda W. Penninx,9 Rick Jansen,9 Kerrie McAloney,4 Jacqueline M. Vink,1 Jaakko Kaprio,10 Robert Plomin,11 Tim D. Spector,12 Patrik K. Magnusson,13 Bruno Reversade,14 R. Alan Harris,15 Kjersti Aagaard,15 Ragnar P. Kristjansson,2 Isleifur Olafsson,16 Gudmundur Ingi Eyjolfsson,17 Olof Sigurdardottir,18 William G. Iacono,5 Cornelis B. Lambalk,19 Grant W. Montgomery,4 Matt McGue,5 Ken K. Ong,20 John R.B Perry,20 Nicholas G. Martin,4 Hreinn Stefánsson,2 Kari Stefánsson,2 and Dorret I. Boomsma,1,21,* 1 Department of Biological Psychology, Vrije U Amsterdamniversiteit, 1081BT Amsterdam, The Netherlands. 2 deCODE Genetics, 101 Reykjavik, Iceland. 3 Institute of Health and Biomedical Innovation, Queensland University of Technology, 4059 Brisbane, Australia. 4 QIMR Berghofer Medical Research Institute, 4029 Brisbane, Australia. 5 Department of Psychology, University of Minnesota, Minneapolis, MN 55455, USA. 6 Queensland Brain Institute, The University of Queensland, 4072 Brisbane, Australia. 7 Avera Institute for Human Genetics, Sioux Falls, SD 57108, USA. 8 Wellcome Trust Centre for Human Genetics, OX3 TBN Oxford, United Kingdom. 9 Department of Psychiatry and EMGO Institute for Health and Care Research, Vrije U Universiteity, Medical Center/GGZ inGeest, 1081HL Amsterdam, The Netherlands. 10 Institute for Molecular Medicine Finland FIMM, University of Helsinki; Department of Public Health, University of Helsinki; National Institute for Health and Welfare, FI-00014 Helsinki, Finland. 1 11 Medical Research Council Social, Genetic and Developmental Psychiatry Centre, Institute of Psychiatry, Psychology & Neuroscience, King's College London, SE5 8AF London, United Kingdom. 12 Department of Twin Research & Genetic Epidemiology, King's College London, SE1 7EH London, United Kingdom. 13 Department of Medical Epidemiology and Biostatistics, Karolinska Institute, SE-171 77 Stockholm, Sweden. 14 Institute of Medical Biology, A*STAR, 138648 Singapore, Singapore. 15 Division of Maternal-Fetal Medicine, Departments of Obstetrics and Gynecology, Molecular and Human Genetics, and Molecular and Cell Biology, Baylor College of Medicine and Texas Children’s Hospital, Houston, TX 77030, USA. 16 Department of Clinical Biochemistry, Landspitali- The National University Hospital of Iceland, 101 Reykjavik, Iceland. 17 Icelandic Medical Center (Laeknasetrid), Laboratory in Mjodd (RAM), Reykjavik, Iceland. 18 Department of Clinical Biochemistry, Akureyri Hospital, 109 Akureyri, Iceland. 19 Department of Obstetrics, Gynecology and Reproductive Medicine, VU University Medical Center, 1007MB Amsterdam, The Netherlands. 20 MRC Epidemiology Unit, University of Cambridge School of Clinical Medicine, Institute of Metabolic Science, Cambridge Biomedical Campus, CB2 0QQ Cambridge, United Kingdom. 21 EMGO+ Institute for Health and Care Research, 1081BT Amsterdam, The Netherlands. 22 These authors contributed equally to this study *Correspondence to h.mbarek@vu.nl (H.M, @drHamdiMbarek), di.boomsma@vu.nl (D.I.B) 2 Abstract Spontaneous dizygotic (DZ) twinning occurs in 1-4% of women, with familial clustering and unknown pathophysiological pathways and genetic origin. DZ twinning may index increased fertility and has distinct health implications for mother and child. We performed the first a GWA study in 1,980 mothers of spontaneous DZ twins and 12,953 controls. Findings were replicated in a large Icelandic cohort and tested for association across a broad range of fertility traits in women. Two SNPs were identified (rs11031006 near FSHB, p = 1.54 × 10-9, and rs17293443 in SMAD3, p = 1.57 × 10-8) and replicated (p = 3 × 10−3 and p = 1.44 × 10−4, respectively). Based on ~90,000 births in Iceland, the relative risk of a mother delivering twins increased by 18% for each copy of the allele rs11031006-G, and 9% for the rs17293443-C. A higher polygenic risk score (PRS) for DZ twinning, calculated based on the results of the DZ twinning GWAS, was significantly associated with DZ twinning in Iceland (p = 0.001). A higher PRS was also associated with having children (p = 0.01), greater lifetime parity (p = 0.03) and earlier age at first child (p = 0.02). The Allele rs11031006-G was associated with higher serum FSH levels, earlier age at menarche, earlier age at first child, higher lifetime parity, lower PCOS risk, and earlier age at menopause. Conversely, the rs17293443-C was associated with later age at last child. We identified the first robust genetic risk variants for DZ twinning: one near FSHB, and a second within SMAD3,—the product of which plays an important role in gonadal responsiveness to FSH. These loci contribute to crucial aspects of reproductive capacity and health. 3 1 Introduction 2 DZ twinning (MIM: 276400) is defined as the concomitant conception and 3 development of two independent zygotes during one pregnancy. Mothers of spontaneous 4 (conceived without assisted reproductive technology) DZ twins have a predisposition to 5 multiple ovulation events due to interference with single dominant follicle selection, a 6 biological mechanism fundamental for the human species.1 DZ twinning is common, with 7 large regional differences from six per 1000 births in Asia to 40 per 1000 births in Africa,2 8 whereas monozygotic twinning occurs around the world at a constant rate of around three to 9 four per 1000 births.3; 4 DZ twinning rates also vary substantially over time. For example, in 10 the US, the observed incidence of twin births increased by a factor of 1.9 between 1971 and 11 2009.5 While a considerable part of the increase is attributable to fertility treatments, with an 12 estimate of 36% of all twins born in the USA in 2011 resulting from assisted reproduction, the 13 majority of twins are still conceived spontaneously. In addition to fertility treatments, 14 increases in maternal age contribute to increases in twinning incidence.3; 6 15 Twinning is associated with increased risks to mother and offspring, including higher 16 risks of stillbirth, neonatal death and premature birth. Compared to singleton children, twins 17 use more hospital resources, especially during the first year of life,7 with hospital costs in the 18 first five years of life being as much as 3·3-fold higher.8 19 Mothers of DZ twins differ from other women in that they are taller, have an increased 20 BMI, are more often overweight and more often smoke before the twin pregnancy.9 Family 21 history, increased parity and gravidity all increase the risk of spontaneous DZ twinning.1; 3 22 Remarkably, twinning rates do not reflect average nutritional status of a population, as 23 established from longitudinal studies in countries that experienced periods of starvation, such 24 as the Dutch hunger winter. Above a specific, yet undetermined, threshold, nutrition seems to 25 be of minimal importance for reproduction in general and also for twinning.10 These 4 26 observations all point to spontaneous twinning being a heritable trait and suggest the potential 27 for polygenic inheritance. In a landmark study of twinning based on data from genealogic 28 records from Utah,11 White and Wyshak established that the genotype of the mother, but not 29 that of the father, affects the frequency of DZ twinning. 30 The underlying physiological mechanism for DZ twinning is the release and 31 fertilization of two or more oocytes. Ovarian folliculogenesis and determination of ovulation 32 quota are controlled both by circulating concentrations of follicle-stimulating hormone (FSH) 33 and by intra-ovarian factors including the two oocyte growth factors, GDF9 (MIM: 601918) 34 and BMP15 (MIM: 300247), as well as their cognate receptors.12 In the common marmoset 35 monkey (subfamily Callitrichinae), DZ twins comprise the predominant litter size, and 36 singletons are rarely, if ever, observed. On DNA sequencing, specific nonsynonymous 37 substitutions were identified in GDF9, BMP15, BMP4 (MIM: 112262), and WFIKKN1 (MIM: 38 608021) as having a role in Callitrichine twinning.13 These genes are among a larger set of 63 39 candidate genes with a potential involvement in regulation of ovulation number and/or control 40 of growth and body size.14 41 Efforts to characterize the genes that contribute to DZ twinning in humans have not 42 been successful. Candidate gene and genome-wide linkage studies failed to uncover common 43 variants associated with DZ twinning,15-19 although one study reported rare variants in GDF9 44 to be associated with DZ twinning.18 45 DZ twinning has been suggested as a measure of human fertility both at the individual 46 and at the population level.20 Spontaneous DZ twinning may be considered as a marker of 47 high fertility, as it reflects the frequency of double ovulation, the probability of coitus within 48 the appropriate time frame with fertilization of both ova, and maintenance of a multiple 49 pregnancy. 5 The aim of this study was to perform the first a genome-wide association study 50 51 (GWAS) in mothers of spontaneous DZ twins to identify relevant genomic regions and test 52 their significance across a broad range of female fertility and reproductive traits including age 53 at menarche, age at natural menopause, age at first and last child, and lifetime parity. Three 54 twin registries, from the Netherlands, Australia and Minnesota (USA) had detailed 55 information on spontaneous twinning in mothers of DZ twins, as well as genotype data. 56 Replication of top hits for twinning was possible in the Icelandic population and for other 57 measures of reproductive ageing in several large-scale population meta-analyses21-24. 58 59 Material and Methods 60 Details of methods and associated references can be found in the supplemental data. 61 Descriptions of Participating Studies 62 Netherlands Twin Register (NTR) 63 The NTR sample consisted of 806 cases and 4,535 controls from the Netherlands Twin 64 Register (2,776 participants) and the Netherlands Study of Depression and anxiety (NESDA; 65 2,565 participants). NTR participants were ascertained by the presence of liveborn twins or 66 triplets in the family and consist of multiples, their parents, siblings and spouses. Twins were 67 born in all strata of society and NTR represents a general sample from the Dutch population25- 68 28 . NESDA is a longitudinal study focusing on the course and consequences of depression and 69 anxiety disorders. Subjects for NESDA were recruited from the general population, mental 70 health organizations and general practices. The sample includes subjects selected for 71 depression and anxiety, as well as healthy controls29. Zygosity of twins was confirmed by 72 DNA genotyping. Data on mode of pregnancy were available from several data collection 73 waves including surveys sent out to mothers of twins, a survey to parents upon registration of 74 young twins, and telephone interview as part of a project on DZ twinning9. The comparison of 6 75 the survey data with the hospital records showed that mothers can accurately report on the 76 mode of conception of their twins30. Participants were excluded if they reported the use of 77 assisted reproductive technology at one or more occasions. In case no reports on mode of 78 pregnancy were available, data were excluded unless the twins were born prior to 1985. 79 QIMR Berghofer Medical Research Institute (QIMR) 80 The sample used in this analysis consisted of 606 cases and 6,656 controls. The individuals 81 were drawn from families containing (any type of) twins recruited for prior studies, either 82 from around Australia from the Australian Twin Registry (ATR) (generally twins born before 83 1971) or from south-eastern Queeensland (generally twins born after 1980). Study recruitment 84 was predominantly population based (any family where the twins were willing to participate) 85 with no screening performed on reproductive phenotypes apart from selecting families with 86 twins. Studies for the older cohort typically were focused on personality traits but not selected 87 for them. Zygosity of twins was reported at time of recruitment; during phenotyping studies; 88 and tested by genotyping (on SNP arrays or Sequenom assays). New phenotyping excluded 89 mothers (‘cases’) from Queensland cohort, who used Assisted Reproductive Technology 90 (ART–typically IVF or hormone treatment) to become pregnant. Screening questions asked of 91 the mother during clinical sessions for phenotyping, were the primary basis for excluding 92 ART cases. A smaller subset of mothers was contacted specifically to establish this 93 information where it was not otherwise available. For the older (ATR) cohort, at the time of 94 the twins’ birth, IVF was not yet in clinical use and other ART was rare. 95 Minnesota Center for Twin and Family Research (MCTFR) 96 All subjects in this sample were independently ascertained through vital records of the State 97 of Minnesota in an effort to construct a population-based twin registry31; 32. The sample for 98 the current study consisted of 568 mothers of DZ twins and 1,862 controls who were the 99 parents of MZ twins from 1,062 families, including 800 complete parental pairs, 203 mothers 7 100 and 59 fathers. Most of the twins were born in the 70s or early 80s, when even though fertility 101 treatment was available in the US, it was expensive and few had access to it. Genotyping was 102 population-based and independent of phenotypes other than twinning. About 92% of the 103 registry, and 100% of both case and control samples, are of primarily European ancestry. 104 Iceland (deCODE) 105 Mothers of twins or other multiples (“twins”) were selected from among those taking part in 106 deCODE’s genetic studies based on a nationwide genealogical database. To increase the 107 proportion of these mother of twins who were mothers of dizygotic twins, twins who had been 108 genotyped and shown to be monozygotic were not used to identify mothers. Controls were 109 individuals participating in deCODE’s genetic research from which both mothers of twins and 110 the mothers’ first-degree relatives had been removed. For the prediction of twinning using 111 polygenic risk scores, mothers having opposite sex or verified dizygotic twins were compared 112 with mothers who did not have twins. 113 114 Study Design, Genome-Wide Association Study and Replication 115 We established the Twinning GWAS Consortium (TGC) to characterize the genetic basis of 116 DZ twinning in humans and performed the first a genome-wide association study (GWAS) 117 for DZ twinning utilizing data from 1,980 mothers of DZ twins (MODZT) and 12,953 118 controls from The Netherlands, US and Australian European ancestry cohorts. Sample sizes 119 and study characteristics are described in table S1. All cases underwent screening to exclude 120 mothers who received assisted reproductive techniques. Controls were screened to exclude 121 pedigrees containing DZ twins. In the replication stage, significant findings from the meta- 122 analysis were tested in large Icelandic cohort of 3,597 mothers of twins and 297,348 controls. 123 All participants provided written informed consent, including consent for genotyping and 8 124 analysis, and were recruited according to the protocols approved by the institution review 125 board of each institution. 126 127 Fertility Traits Measures 128 Fertility measures (having children, number of children, age at first and last child and average 129 birth interval) were defined in females born before 1970 based on the Icelandic genealogy. 130 Data for age at menarche, age at menopause and polycystic ovary syndrome were derived 131 from previously published GWAS consortia21-24. Sample sizes and study characteristics are 132 described in table 2. 133 134 Genotyping, Quality Control and Imputation 135 Each participating cohort performed participant-level genotyping of single nucleotide 136 polymorphisms (SNPs), that included standard quality-control measures for genotyping and 137 imputation (Table S21 and Supplemental Methods). All imputations were performed using the 138 1000 Genomes Project March 2012 release as the reference panel. 139 Statistical analysis 140 GWA analysis was performed in each cohort using logistic regression under an additive 141 genetic model with adjustment for principal components of genetic ancestry and specific 142 covariates for each study. GWA results for each SNP (odds ratio [OR] and 95% confidence 143 interval [CI]) were combined across cohorts using fixed-effect meta-analysis with inverse 144 variance weighting. The most significant SNP at each genome-wide significant (p<5×10-8) 145 locus was tested in the replication stage and replicating SNPs were examined for association 146 with FSH levels and other fertility traits. 147 Association Tests 9 148 Each genome-wide association analysis from the three cohorts was conducted using logistic 149 regression under an additive genetic model with adjustment for principal components of 150 genetic ancestry. Because the GWAS data include family members we added the –family 151 option in the analysis, which takes the familial structure of the data into account using a 152 sandwich estimator33. Imputed SNPs were analyzed using PLINK software34 and genotype 153 imputation uncertainty was accounted for by using allelic dosage in PLINK. 154 Meta-analysis was performed using the fixed-effects inverse variance method based on the 155 regression β estimates and standard errors from each study implemented in METAL35. The 156 presence of heterogeneity between cohorts for the effect sizes of risk alleles was investigated 157 using the Cochran’s Q-test as implemented in METAL. To determine whether the genome- 158 wide significant signal at each locus with low LD in the same chromosomal region (defined 159 as r2 < 0.05 in a 750-kb region) could be accounted for by a single SNP, we carried out 160 conditional analysis. Each cohort performed a genome-wide analysis for MODZT using 161 logistic regression adjusting for the top signal at each of the three associated regions to 162 determine whether potential second signals remained significant even after adjusting for these 163 variants. Results from each individual study were meta-analyzed to determine whether these 164 potential second signals were truly independent (that is, if p < 5 × 10-8). In Iceland, being a 165 mother of twins was tested for association with the top alleles using logistic regression and 166 including age, age-squared and county of birth as covariates as described previously36. 167 Associations between FSH level and genotype was assessed using linear regression as 168 described previously36. To study the association of fertility measures and SNP genotype, we 169 used logistic regression (having children), Poisson log-linear regression (number of children), 170 or linear regression (age at first child, age at last child, and average birth interval). In these 171 analyses, birth cohort (as a factor for each five year interval), county of birth and six principal 172 components were included as covariates. Relatedness was controlled for by using genomic 10 173 control in all Icelandic association analyses37. The combined p value of the meta-analysis and 174 the replication study was calculated using Fisher’s combined probability test 38. 175 176 FSH Serum Level Measure 177 Serum levels of follicle-stimulating hormone (FSH) were measured in 2,411 men (1,275 178 genotyped persons; 1,136 close relatives of genotyped individuals) and 15,586 women (9,738 179 genotyped persons; 5,848 close relatives of genotyped individuals) referred to three clinics in 180 Iceland. FSH testing was undertaken primarily to investigate possible gonad impairment. 181 Hormone levels were measured by electrochemiluminescence immunoassay, using reagents 182 and analytical instruments from Roche Diagnostics GmbH, according to the manufacturer’s 183 instructions. 184 185 186 In silico Functional Annotation We used a number of publicly available bioinformatics tools and datasets to identify putative 187 functional effects of the top associated SNPs at each locus, including: Combined Annotation 188 Dependent Depletion (CADD)39, HaploReg2 v440 and Variant effect predictor (VEP)41 . 189 190 Gene Based Test 191 Results of the MODZT meta-analysis were used to perform a gene-based test of association 192 for the 63 candidate genes from Harris et al. using the Knowledge-based mining system for 193 Genome-wide Genetic studies (KGG) software Version 3.542; 43. This approach uses an 194 extended Simes test that integrates prior functional information and the meta-analysis 195 association results when combining the SNP p values within a gene in order to obtain an 196 overall association p value for each entire gene. As we tested for genetic association for 63 197 genes, the significance level was set at 7.93 × 10-4 (Bonferroni correction; 0.05/63). 11 198 199 Polygenic Risk Scores (PRS) 200 Polygenic risk scores were calculated based on the results of the MODZT meta-analysis. Only 201 markers having info > 0.9 in all groups and MAF > 0.01 were included. To obtain effect sizes 202 taking LD into account, the LDpred method developed by Vilhjálmsson and colleagues was 203 used44. As suggested by Vilhjálmsson et al, we calculated multiple sets of LD-modified effect 204 sizes based on a grid of values for the fraction of causal markers (α = 0.0001, 0.0003, 0.001, 205 0.003, 0.01, 0.03, 0.1, 0.3, 1). The resulting scores were then tested in a validation data set of 206 Icelandic mothers of dizygotic twins. The score producing the most significant result in the 207 validation data set were subsequently used to test for association with five fertility-related 208 traits (“has children”, “number of children”, “age at first child”, “age at last child” and 209 “average birth interval”). 210 211 Results 212 Genome-wide association results, replication and gene based analyses 213 The overall GWAS meta-analysis genomic control statistic (λ) was 1.01, indicating no 214 appreciable inflation due to population structure. The quantile-quantile (Q-Q) plot of genome- 215 wide p values showed a strong deviation from the null hypothesis of no association (Figure 216 S1). The results are represented in the Manhattan plot (Figure 1). Three chromosomal regions 217 contained genome-wide significant SNPs (p < 5 × 10-8). Twenty-two SNPs on chromosome 218 11p13 showed genome-wide significant associations. Of these, the strongest signal 219 (rs11031006, p = 1.54 × 10-9 and OR 1.41, 95% CI 1.29-1.53) lies in the region 5’ of the 220 transcription start site of FSHB (encoding the FSH, beta polypeptide, [MIM: 136530]). A 221 second locus was represented by five SNPs on 1q42.13 covering an intergenic region flanked 222 by PGBD5 (MIM: 616791) and COG2 ([MIM: 606974]), (p = 1.23 × 10-8, OR 1.43, 95% CI 12 223 1.31-1.55 for strongest signal, rs12064669). A third locus on 15q22.33, mapped to the first 224 intron of SMAD3 (SMAD family member 3, [MIM: 603109]) (rs17293443, p = 1.57 × 10-8 225 and OR 1.27, 95% CI 1.19-1.35). No significant heterogeneity in SNP effects was observed 226 across cohorts for the top SNPs (p > 0.1, Cochran’s Q test) ; Tables 1 3 and S23). The 227 regional association plots for these loci are shown in figure S2. After conditioning on the top 228 SNPs at each locus, no secondary signals were observed (all p > 0.05; Tables S43, S54 and 229 S65). We sought validation of these three top signals in an independent replication study from 230 Iceland (deCODE) totaling 3,597 mothers of twins and 297,348 controls. The FSHB 231 (rs11031006, p = 3 × 10−3, OR 1.14, 95% CI 1.06-1.22) and SMAD3 (rs17293443, p = 1.44 × 232 10−4, OR 1.15, 95% CI 1.07-1.23) loci replicated, but rs12064669 (p=0.88) was not confirmed 233 (Table 13 and Figure S32). We also investigated whether any of the proposed 63 candidate 234 genes14 was associated with human DZ twinning. In a gene-based test, five genes 235 demonstrated a nominally significant association (p < 0.05; BMPR1A [MIM: 601299], 236 BMPR1B [MIM: 603248], IGF1 [MIM: 147440], FSHB and FSHR [MIM: 136435]). After 237 correction for multiple testing, only FSHB remained significant (Table S76). 238 239 Functional in silico annotation of associated variants 240 We explored plausible functional effects of our associated variants using Combined 241 Annotation Dependent Depletion (CADD).39 The FSHB SNP rs11031006 had a high Phred 242 scaled C-score (22.4), indicating that it is among the top 1% of SNPs in the human genome 243 most likely to have a functional effect. The Phred score for the SMAD3 SNP rs17293443 was 244 only 2.71 indicating that it is among the bottom 50% of SNPs in the human genome likely to 245 have a functional effect (Table S107). Examination of individual constituents of the CADD 246 scores showed particularly high conservation-based scores for rs11031006 (Figure S63). To 247 further investigate possible functional effects we examined data from the ENCODE project. 45 Field Code Changed 13 Field Code Changed 248 The FSHB SNP rs11031006 alters the sequence of 11 protein-binding motifs including that of 249 the Estrogen receptor alpha, indicating a possible effect on hormonal feedback inhibition 250 (Table S118). Variant Effect Predictor (VEP)41 identified rs17293443 as a regulatory region 251 variant within a promoter flanking region (ENSR00000410126) (Figure S74). SNP 252 rs17293443 was contained in a DNase I hypersensitive site suggesting open chromatin, 253 although it did not alter any of the transcription factor binding sites present in the promoter 254 flanking region. Together, these data indicate that rs11031006 and rs17293443 may have 255 direct functional roles. 256 SNPs associated with FSH serum level 257 We analyzed serum FSH measurements from 17,997 genotyped Icelanders and their close 258 relatives. SNP rs11031006 was significantly associated with higher serum FSH levels (p = 2.3 259 × 10-10), with each G allele conferring an increase in FSH level of about 0.11 SD units (Figure 260 S54). Notably, the allele (rs11031006-G) conferring the strongest association with FSH levels 261 in the FSH GWAS is the same allele that conferred the greatest chance of having DZ twins in 262 our MODZT GWAS. No association was seen between the SMAD3 signal and serum FSH 263 levels (p = 0.30). 264 265 Relative risk of twin birth and polygenic risk score prediction of DZ twinning and 266 fertility measures 267 We estimated that the relative risk of a twin birth, based on approximately 90,000 births in 268 Iceland between 1950 and 1991, was increased by 18% for each maternal rs11031006-G 269 allele and by 9% for each rs17293443-C allele (Table S8). A higher polygenic risk score 270 (PRS) for DZ twinning, calculated based on the results of the DZ twinning GWAS, was 271 significantly associated with DZ twinning in our independent Icelandic cohort (p = 0.001) 272 (Table S7). A higher PRS was also associated with a higher likelihood of having children (p = 14 273 0.01), higher lifetime number of children (p = 0.03), and an earlier age at first child (p = 0.02) 274 (Table S98). A re-calculated PRS, excluding the 1 Mb regions surrounding the two replicated 275 variants, remained associated with DZ twinning (p = 0.02) and with the likelihood of having 276 children (p = 0.03). These results reflect the polygenic contribution to the susceptibility to DZ 277 twinning and its association with greater reproductive ability. 278 279 SNPs associated with female reproduction traits 280 Table 24 reports on the two loci robustly implicated in DZ twinning and other reproductive 281 traits in women. Consistent with its effects on higher circulating FSH levels, the rs11031006- 282 G allele is also associated with earlier age at menarche,22; 46 earlier age at first child and 283 higher total lifetime number of children, lower risk of polycystic ovary syndrome (PCOS),23 284 and earlier age at natural menopause.24; 47 Also, the DZ twinning SNP rs11031006 is 285 correlated with a variant (rs10835638, FSHB-211G>T) located in the promoter of FSHB (r2 = 286 0.62) that is associated with timing of breast development in girls.48 In contrast, the 287 rs17293443-C allele in SMAD3 was associated only with a later age at last child (Figure S65). 288 Functional in silico annotation of associated variants 289 We explored plausible functional effects of our associated variants using Combined 290 Annotation Dependent Depletion (CADD).38 The FSHB SNP rs11031006 had a high Phred 291 scaled C-score (22.4), indicating that it is among the top 1% of SNPs in the human genome 292 most likely to have a functional effect. The Phred score for the SMAD3 SNP rs17293443 was 293 only 2.71 indicating that it is among the bottom 50% of SNPs in the human genome likely to 294 have a functional effect (Table S10). Examination of individual constituents of the CADD 295 scores showed particularly high conservation-based scores for rs11031006 (Figure S6). To 296 further investigate possible functional effects we examined data from the ENCODE project.44 297 The FSHB SNP rs11031006 alters the sequence of 11 protein-binding motifs including that of Field Code Changed 15 Field Code Changed 298 the Estrogen receptor alpha, indicating a possible effect on hormonal feedback inhibition 299 (Table S11). Variant Effect Predictor (VEP)40 identified rs17293443 as a regulatory region 300 variant within a promoter flanking region (ENSR00000410126) (Figure S7). SNP rs17293443 301 was contained in a DNase I hypersensitive site suggesting open chromatin, although it did not 302 alter any of the transcription factor binding sites present in the promoter flanking region. 303 Together, these data indicate that rs11031006 and rs17293443 may have direct functional 304 roles. 305 Discussion 306 Here we report the first compelling evidence that sequence variation at the FSHB and 307 SMAD3 loci increases the odds of DZ twinning. A number of studies in mothers of DZ twins, 308 but not all,49 have found higher FSH levels responsible for multiple follicle growth.50 The 309 associations of FSHB rs11031006-G with earlier ages at breast development, menarche, 310 menopause, first child, and higher lifetime parity indicates that this locus plays an important 311 role in multiple reproductive aspects. Female carriers of rs11031006-G likely may have a 312 more advanced depletion of the ovarian follicular pool and hence would have an increased 313 risk of premature ovarian failure (POF, [MIM: 612964]). Indeed advanced ovarian aging is a 314 recognized feature of familial DZ twinning, with reported lower levels of anti-Mullerian 315 Hormone (AMH) a marker of lower ovarian primordial follicular reserve.51 The rs17293443- 316 C allele in SMAD3 also increases chances of DZ twinning, but this effect appears independent 317 of circulating FSH levels. 318 Of the 63 suggested candidate genes for twinning, only FSHB was associated in gene- 319 based tests after correcting for multiple testing. Emerging data in human and non-human 320 primates suggest that mechanisms underlying multiple ovulation may differ among species, 321 explaining why in some species twinning is accompanied by unique evolutionary adaptations 322 enabling offspring’s survival, such as the marmoset diminutive fetal size in a simplex uterus, 16 323 and subsequent alloparenting.14 Conversely, in humans, multiple gestations remain an 324 independent risk factor for preterm birth, pregnancy loss, and fetal growth restriction. Thus 325 loci common to species and strains with higher rates of DZ twinning will not necessarily be 326 shared, as the rate, complications, and future reproductive fitness of those twin gestations 327 differ. 328 The third genome-wide hit from the discovery was not replicated in Iceland. SNP 329 rs120644669 is located 149kb from the angiotensinogen (AGT, [MIM: 106150]), which 330 influences ovulatory capacity in mice,52 and 89kb from component of oligomeric golgi 331 complex 2 (COG2) affecting protein glycosilation, which regulates the biological activity of 332 the pituitary gonadotrophins.53 333 Genetic variants near FSHB (rs11031005 and rs11031002, which are highly correlated 334 with rs11031006) are reportedly associated not only with higher serum FSH, but also lower 335 LH levels.54 This agrees with the known response of the gonadotrophic cell to a high 336 frequency hypothalamic pulsatile GnRH signal that reciprocally controls secretion of both 337 hormones.55 Furthermore, under normal conditions suppression of FSH in the early follicular 338 phase and higher LH levels in late phase typically favor mono-ovulation in the human.56 In a 339 study aiming to identify genetic predictors for IVF success or IVF-controlled ovarian 340 stimulation (COS), rs611246 located in FSHB (r2 = 0.3 with rs11031006), was reported 341 significantly associated with measured early follicular phase FSH values and also with the 342 probability of clinical pregnancy, suggesting that these genetic variants are potential predictor 343 candidates that could be considered in clinical ovarian reserve and function assessment in 344 assisted reproduction.57 345 A recent linkage study in cattle reported only one strong signal (p < 1 × 10-28) for 346 ovulation rate in a region spanning SMAD3, SMAD6 (MIM: 602931) and IQCH (MIM: 347 612523).58 SMAD3 encodes one of a family of proteins that function as signal transducers and 17 348 transcriptional modulators that mediate multiple signaling pathways. Observations in mice 349 have established an essential role for Smad3 in mediating TGFbeta and activin signals in the 350 ovarian granulosa cell and also in the pituitary59 to maintain a favorable environment for 351 oocyte maturation.60 Smad3 is strongly expressed in the human ovary, where it promotes 352 granulosa cell proliferation and steroidogenesis possibly by upregulating gonadotrophin 353 receptor signaling pathways.61 Thus, sequence variation in SMAD3 may increase chance of 354 DZ twinning by increasing responsiveness to FSH. Understanding the role of SMAD3 will 355 offer novel opportunities to optimize responsiveness and minimize risk among assisted 356 reproduction technology (ART) recipients for example through adjustment of hormonal 357 stimulation and thus contributes to prevention of life-threatening ovarian hyper stimulation 358 syndrome in hyper responding female carriers of the rs17293443-C allele and conversely 359 prevention of a poor response in patients individuals with allelic variants that lead to a poor 360 response. 361 The potential applications of this work in reproductive medicine are multifold. Firstly, 362 it reveals a well-defined set of loci for DZ twinning in the general population. Secondly, 363 twinning is associated with common perinatal morbidities such as preterm birth, discordant 364 twin growth, latter prenatal asymmetric intrauterine growth restriction, and placental 365 abruption. Multiple gestations are also related to a higher prevalence of maternal morbidities 366 such as preeclampsia, postpartum hemorrhage, and ensuing complications. By understanding 367 the genetic basis of DZ twinning, we concomitantly identify loci conferring susceptibility (or 368 conversely resistance) to these prevalent perinatal comorbidities. The involvement of SMAD3 369 is of relevance in the light of a possible phenotype in male DZ twins namely the repeatedly 370 reported more frequent occurrence of testicular seminoma, a gonadal tumor in which the 371 mitogenic cyclin D2 is overexpressed.62; 63 In animal experimental studies knockout of 372 SMAD3 significantly attenuates cyclin D2 and tumor proliferation.64 This indicates to the 18 373 likely clinical potential of SMAD3 stretching beyond that of fertility physiology and 374 treatment. 375 This study The sequence variants found to influence spontaneous DZ twinning and 376 their relationship with other fertility-related measures, provides important insights into 377 ovarian functioning and the control of natural multiple follicle growth and reproductive aging. 378 This has important implications for fertility, including improved outcome prediction and 379 novel avenues of fertility treatment. Other strengths include the rigorous phenotype inclusion 380 of mothers of DZ twins with a documented history of spontaneous DZ twinning. It is worth 381 mentioning that analyses without excluding mothers who conceived their twins with hormone 382 induction of multiple ovulation or other ART, did not yield to any genome-wide significant 383 results (results not shown). We thus recommend twins registries collect data on mode of 384 pregnancy of twins. There are also some limitations. The current effort focused on unraveling 385 the genetic basis of DZ twinning in European-ancestry populations only. However, the 386 highest incidence of DZ twins was reported in the Nigerian population and some other 387 countries from Africa. Next steps in unraveling the genetic cause of DZ twinning need to 388 include mothers of DZ twins originating from these regions, and if possible also from regions 389 where DZ twinning is a rare trait, such as Japan2. 390 In summary, we identified FSHB and the novel SMAD3 locus as maternal 391 susceptibility loci for DZ twinning. These loci are also significantly and specifically 392 associated with several other aspects of reproductive capacity and health. 393 394 395 396 19 397 398 399 Supplemental Data Supplemental data include seven six figures and 11eight tables Acknowledgments 400 A list of support provided to individual studies appears in supplemental data. Support for the 401 Netherlands Twin Register was obtained from the Netherlands Organization for Scientific 402 Research (NWO) and The Netherlands Organisation for Health Research and Development 403 (ZonMW) grants, 904-61-193,480-04-004, 400-05-717, Addiction-31160008, 911-09-032, 404 Spinozapremie 56-464-14192, Biobanking and Biomolecular Resources Research 405 Infrastructure (BBMRI –NL, 184.021.007); the European Research Council (ERC-230374 406 and ERC-284167); Rutgers University Cell and DNA Repository (NIMH U24 MH068457- 407 06), the Avera Institute, Sioux Falls, South Dakota (USA) and the National Institutes of 408 Health (NIH R01 HD042157-01A1). Part of the genotyping was funded by the Genetic 409 Association Information Network (GAIN) of the Foundation for the National Institutes of 410 Health and Grand Opportunity grants 1RC2 MH089951). We acknowledge support from VU 411 Amsterdam and the Institute for Health and Care Research (EMGO+). The Berghofer Medical 412 Research Institute (QIMR) study was supported by grants from the National Health and 413 Medical Research Council (NHMRC) of Australia (241944, 339462, 389927, 389875, 414 389891, 389892, 389938, 443036, 442915, 442981, 496610, 496739, 552485, 552498, 415 1050208, 1075175). Dale R. Nyholt was supported by the Australian Research Council 416 (ARC) Future Fellowship (FT0991022), NHMRC Research Fellowship (APP0613674) 417 Schemes and by the Visiting Professors Programme (VPP) of the Royal Netherlands 20 418 Academy of Arts and Sciences (KNAW). Allan F. McRae was supported by an NRMRC 419 Career Development Fellowship (APP1083656). Grant W. Montgomery was supported by 420 NIH grant (HD042157, a collaborative study of the genetics of DZ twinning) and NHMRC 421 Fellowship (GNT1078399). The Minnesota Center for Twin and Family Research (MCTFR) 422 was supported in part by USPHS Grants from the National Institute on Alcohol Abuse and 423 Alcoholism (AA09367 and AA11886), and the National Institute on Drug Abuse (DA05147, 424 DA13240, and DA024417). 425 We would like to thank also 23andMe and 23andMe's consented research participants for 426 contributing data on age at menarche for the FSHB gene locus and the Twinning Gwas 427 Consortium (TGC) co-authros from: Finland (Anu Loukola, Juho Wedenoja, Emmi Tikkanen, 428 Beenish Qaiser), Sweden (Nancy Pedersen, Andrea Ganna), United kingdom King's College 429 London (Department of Twin Research & Genetic Epidemiology: Pirro Hysi, Massimo 430 Mangino), Institute of Psychiatry, Psychology & Neuroscience, Medical Research Council 431 Social, Genetic and Developmental Psychiatry Centre (Eva Krapohl, Andrew McMillan). S.S, R.P.K, H.S and K.S are employees of deCODE Genetics/Amgen. The other authors declare no competing financial interests. Web Resources 1000GenomesProject, ftp://ftp-trace.ncbi.nih.gov/1000genomes/ ftp/release/20110521/ CADD, http://cadd.gs.washington.edu/ HaploReg, http://www.broadinstitute.org/mammals/haploreg/haploreg.php Variant Effect Predictor (VEP), http://www.ensembl.org/info/docs/tools/vep/index.html Knowledge-based mining system for Genome-wide Genetic studies (KGG), http://grass.cgs.hku.hk/limx/kgg/ ENCODE, https://genome.ucsc.edu/ENCODE/ OMIM, http://www.omim.org/ 21 LDpred, https://bitbucket.org/bjarni_vilhjalmsson/ldpred MACH, http://www.sph.umich.edu/csg/abecasis/MACH/ Minimac, http://genome.sph.umich.edu/wiki/Minimac Beagle, https://faculty.washington.edu/browning/beagle/b3.html PLINK, http://pngu.mgh.harvard.edu/~purcell/plink/ LocusZoom, http://csg.sph.umich.edu/locuszoom/ References 432 1. Bulmer, M.G. (1970). The biology of twinning in man.(Oxford,: Clarendon). 433 2. Hall, J.G. (2003). Twinning. Lancet 362, 735-743. 434 435 3. Hoekstra, C., Zhao, Z.Z., Lambalk, C.B., Willemsen, G., Martin, N.G., Boomsma, D.I., and Montgomery, G.W. (2008). Dizygotic twinning. Hum Reprod Update 14, 37-47. 436 437 4. Smits, J., and Monden, C. (2011). Twinning across the Developing World. PLoS One 6, e25239. 438 439 440 5. Kulkarni, A.D., Jamieson, D.J., Jones, H.W., Jr., Kissin, D.M., Gallo, M.F., Macaluso, M., and Adashi, E.Y. (2013). Fertility treatments and multiple births in the United States. N Engl J Med 369, 2218-2225. 441 442 6. Fauser, B.C., Devroey, P., and Macklon, N.S. (2005). Multiple birth resulting from ovarian stimulation for subfertility treatment. Lancet 365, 1807-1816. 443 444 445 446 7. Chambers, G.M., Hoang, V.P., Lee, E., Hansen, M., Sullivan, E.A., Bower, C., and Chapman, M. (2014). Hospital costs of multiple-birth and singleton-birth children during the first 5 years of life and the role of assisted reproductive technology. JAMA Pediatr 168, 1045-1053. 447 448 449 450 8. van Heesch, M.M., Evers, J.L., van der Hoeven, M.A., Dumoulin, J.C., van Beijsterveldt, C.E., Bonsel, G.J., Dykgraaf, R.H., van Goudoever, J.B., Koopman-Esseboom, C., Nelen, W.L., et al. (2015). Hospital costs during the first 5 years of life for multiples compared with singletons born after IVF or ICSI. Hum Reprod 30, 1481-1490. 451 452 453 9. Hoekstra, C., Willemsen, G., van Beijsterveldt, C.E., Lambalk, C.B., Montgomery, G.W., and Boomsma, D.I. (2010). Body composition, smoking, and spontaneous dizygotic twinning. Fertil Steril 93, 885-893. 454 455 456 10. Eriksson, A.W., Bressers, W.M., Kostense, P.J., Pitkanen, K.J., Mielke, J.H., Jorde, L.B., Tas, R.F., and Fellman, J.O. (1988). Twinning rate in Scandinavia, Germany and The Netherlands during years of privation. Acta Genet Med Gemellol (Roma) 37, 277-297. 457 458 11. White, C., and Wyshak, G. (1964). Inheritance in Human Dizygotic Twinning. N Engl J Med 271, 1003-1005. 459 460 12. Moore, R.K., Erickson, G.F., and Shimasaki, S. (2004). Are BMP-15 and GDF-9 primary determinants of ovulation quota in mammals? Trends Endocrinol Metab 15, 356-361. 22 461 462 13. (2014). The common marmoset genome provides insight into primate biology and evolution. Nat Genet 46, 850-857. 463 464 465 14. Harris, R.A., Tardif, S.D., Vinar, T., Wildman, D.E., Rutherford, J.N., Rogers, J., Worley, K.C., and Aagaard, K.M. (2014). Evolutionary genetics and implications of small size and twinning in callitrichine primates. Proc Natl Acad Sci U S A 111, 1467-1472. 466 467 15. Lambalk, C.B. (2001). Is there a role for follicle-stimulating-hormone receptor in familial dizygotic twinning? Lancet 357, 735-736. 468 469 470 16. Duffy, D., Montgomery, G., Treloar, S., Birley, A., Kirk, K., Boomsma, D., Beem, L., de Geus, E., Slagboom, E., Knighton, J., et al. (2001). IBD sharing around the PPARG locus is not increased in dizygotic twins or their mothers. Nat Genet 28, 315. 471 472 473 474 17. Painter, J.N., Willemsen, G., Nyholt, D., Hoekstra, C., Duffy, D.L., Henders, A.K., Wallace, L., Healey, S., Cannon-Albright, L.A., Skolnick, M., et al. (2010). A genome wide linkage scan for dizygotic twinning in 525 families of mothers of dizygotic twins. Hum Reprod 25, 1569-1580. 475 476 477 478 18. Palmer, J.S., Zhao, Z.Z., Hoekstra, C., Hayward, N.K., Webb, P.M., Whiteman, D.C., Martin, N.G., Boomsma, D.I., Duffy, D.L., and Montgomery, G.W. (2006). Novel variants in growth differentiation factor 9 in mothers of dizygotic twins. J Clin Endocrinol Metab 91, 4713-4716. 479 480 481 482 19. Luong, H.T., Chaplin, J., McRae, A.F., Medland, S.E., Willemsen, G., Nyholt, D.R., Henders, A.K., Hoekstra, C., Duffy, D.L., Martin, N.G., et al. (2011). Variation in BMPR1B, TGFRB1 and BMPR2 and control of dizygotic twinning. Twin Res Hum Genet 14, 408-416. 483 484 20. Tong, S., Caddy, D., and Short, R.V. (1997). Use of dizygotic to monozygotic twinning ratio as a measure of fertility. Lancet 349, 843-845. 485 486 487 21. Day, F.R., Bulik-Sullivan, B., Hinds, D.A., Finucane, H.K., Murabito, J.M., Tung, J.Y., Ong, K.K., and Perry, J.R. (2015). Shared genetic aetiology of puberty timing between sexes and with health-related outcomes. Nat Commun 6, 8842. 488 489 490 22. Perry, J.R., Day, F., Elks, C.E., Sulem, P., Thompson, D.J., Ferreira, T., He, C., Chasman, D.I., Esko, T., Thorleifsson, G., et al. (2014). Parent-of-origin-specific allelic associations among 106 genomic loci for age at menarche. Nature 514, 92-97. 491 492 493 494 23. Day, F.R., Hinds, D.A., Tung, J.Y., Stolk, L., Styrkarsdottir, U., Saxena, R., Bjonnes, A., Broer, L., Dunger, D.B., Halldorsson, B.V., et al. (2015). Causal mechanisms and balancing selection inferred from genetic associations with polycystic ovary syndrome. Nat Commun 6, 8464. 495 496 497 498 24. Day, F.R., Ruth, K.S., Thompson, D.J., Lunetta, K.L., Pervjakova, N., Chasman, D.I., Stolk, L., Finucane, H.K., Sulem, P., Bulik-Sullivan, B., et al. (2015). Large-scale genomic analyses link reproductive aging to hypothalamic signaling, breast cancer susceptibility and BRCA1-mediated DNA repair. Nat Genet 47, 1294-1303. 499 500 501 502 25. Boomsma, D.I., de Geus, E.J., Vink, J.M., Stubbe, J.H., Distel, M.A., Hottenga, J.J., Posthuma, D., van Beijsterveldt, T.C., Hudziak, J.J., Bartels, M., et al. (2006). Netherlands Twin Register: from twins to twin families. Twin research and human genetics : the official journal of the International Society for Twin Studies 9, 849-857. 503 504 26. Willemsen, G., de Geus, E.J., Bartels, M., van Beijsterveldt, C.E., Brooks, A.I., Estourgievan Burk, G.F., Fugman, D.A., Hoekstra, C., Hottenga, J.J., Kluft, K., et al. (2010). 23 505 506 The Netherlands Twin Register biobank: a resource for genetic epidemiological studies. Twin Res Hum Genet 13, 231-245. 507 508 509 510 27. van Beijsterveldt, C.E., Groen-Blokhuis, M., Hottenga, J.J., Franic, S., Hudziak, J.J., Lamb, D., Huppertz, C., de Zeeuw, E., Nivard, M., Schutte, N., et al. (2013). The Young Netherlands Twin Register (YNTR): longitudinal twin and family studies in over 70,000 children. Twin Res Hum Genet 16, 252-267. 511 512 513 514 28. Willemsen, G., Vink, J.M., Abdellaoui, A., den Braber, A., van Beek, J.H., Draisma, H.H., van Dongen, J., van 't Ent, D., Geels, L.M., van Lien, R., et al. (2013). The Adult Netherlands Twin Register: twenty-five years of survey and biological data collection. Twin Res Hum Genet 16, 271-281. 515 516 517 518 29. Penninx, B.W., Beekman, A.T., Smit, J.H., Zitman, F.G., Nolen, W.A., Spinhoven, P., Cuijpers, P., De Jong, P.J., Van Marwijk, H.W., Assendelft, W.J., et al. (2008). The Netherlands Study of Depression and Anxiety (NESDA): rationale, objectives and methods. International journal of methods in psychiatric research 17, 121-140. 519 520 521 522 30. van Beijsterveldt, C.E., Hoekstra, C., Schats, R., Montgomery, G.W., Willemsen, G., and Boomsma, D.I. (2008). Mode of conception of twin pregnancies: willingness to reply to survey items and comparison of survey data to hospital records. Twin Res Hum Genet 11, 349-351. 523 524 525 526 31. Miller, M.B., Basu, S., Cunningham, J., Eskin, E., Malone, S.M., Oetting, W.S., Schork, N., Sul, J.H., Iacono, W.G., and McGue, M. (2012). The Minnesota Center for Twin and Family Research genome-wide association study. Twin Res Hum Genet 15, 767774. 527 528 529 32. Iacono, W.G., Carlson, S.R., Taylor, J., Elkins, I.J., and McGue, M. (1999). Behavioral disinhibition and the development of substance-use disorders: findings from the Minnesota Twin Family Study. Dev Psychopathol 11, 869-900. 530 531 33. Rogers, W.H. (1993). Regression standard errors in clustered samples. Stata Technical Bulletin 13, 19-23. 532 533 534 535 34. Purcell, S., Neale, B., Todd-Brown, K., Thomas, L., Ferreira, M.A., Bender, D., Maller, J., Sklar, P., de Bakker, P.I., Daly, M.J., et al. (2007). PLINK: a tool set for wholegenome association and population-based linkage analyses. Am J Hum Genet 81, 559575. 536 537 35. Willer, C.J., Li, Y., and Abecasis, G.R. (2010). METAL: fast and efficient meta-analysis of genomewide association scans. Bioinformatics 26, 2190-2191. 538 539 540 541 36. Gudbjartsson, D.F., Helgason, H., Gudjonsson, S.A., Zink, F., Oddson, A., Gylfason, A., Besenbacher, S., Magnusson, G., Halldorsson, B.V., Hjartarson, E., et al. (2015). Large-scale whole-genome sequencing of the Icelandic population. Nat Genet 47, 435444. 542 543 37. Devlin, B., and Roeder, K. (1999). Genomic control for association studies. Biometrics 55, 997-1004. 544 545 38. Fisher, R.A. (1925). Statistical methods for research workers.(Edinburgh, London,: Oliver and Boyd). 546 547 548 39. Kircher, M., Witten, D.M., Jain, P., O'Roak, B.J., Cooper, G.M., and Shendure, J. (2014). A general framework for estimating the relative pathogenicity of human genetic variants. Nat Genet 46, 310-315. 24 549 550 551 40. Ward, L.D., and Kellis, M. (2012). HaploReg: a resource for exploring chromatin states, conservation, and regulatory motif alterations within sets of genetically linked variants. Nucleic Acids Res 40, D930-934. 552 553 554 41. McLaren, W., Pritchard, B., Rios, D., Chen, Y., Flicek, P., and Cunningham, F. (2010). Deriving the consequences of genomic variants with the Ensembl API and SNP Effect Predictor. Bioinformatics 26, 2069-2070. 555 556 557 42. Li, M.X., Kwan, J.S., and Sham, P.C. (2012). HYST: a hybrid set-based test for genomewide association studies, with application to protein-protein interaction-based association analysis. Am J Hum Genet 91, 478-488. 558 559 560 43. Li, M.X., Gui, H.S., Kwan, J.S., and Sham, P.C. (2011). GATES: a rapid and powerful gene-based association test using extended Simes procedure. Am J Hum Genet 88, 283-293. 561 562 563 564 44. Vilhjalmsson, B.J., Yang, J., Finucane, H.K., Gusev, A., Lindstrom, S., Ripke, S., Genovese, G., Loh, P.R., Bhatia, G., Do, R., et al. (2015). Modeling Linkage Disequilibrium Increases Accuracy of Polygenic Risk Scores. Am J Hum Genet 97, 576-592. 565 566 45. (2012). An integrated encyclopedia of DNA elements in the human genome. Nature 489, 57-74. 567 568 569 46. Lunetta, K.L., Day, F.R., Sulem, P., Ruth, K.S., Tung, J.Y., Hinds, D.A., Esko, T., Elks, C.E., Altmaier, E., He, C., et al. (2015). Rare coding variants and X-linked loci associated with age at menarche. Nat Commun 6, 7756. 570 571 572 573 47. Stolk, L., Perry, J.R., Chasman, D.I., He, C., Mangino, M., Sulem, P., Barbalic, M., Broer, L., Byrne, E.M., Ernst, F., et al. (2012). Meta-analyses identify 13 loci associated with age at menopause and highlight DNA repair and immune pathways. Nat Genet 44, 260-268. 574 575 576 577 48. Hagen, C.P., Sorensen, K., Aksglaede, L., Mouritsen, A., Mieritz, M.G., Tinggaard, J., Wohlfart-Veje, C., Petersen, J.H., Main, K.M., Rajpert-De Meyts, E., et al. (2014). Pubertal onset in girls is strongly influenced by genetic variation affecting FSH action. Sci Rep 4, 6412. 578 579 580 49. Gilfillan, C.P., Robertson, D.M., Burger, H.G., Leoni, M.A., Hurley, V.A., and Martin, N.G. (1996). The control of ovulation in mothers of dizygotic twins. J Clin Endocrinol Metab 81, 1557-1562. 581 582 583 584 50. Lambalk, C.B., Boomsma, D.I., De Boer, L., De Koning, C.H., Schoute, E., PoppSnijders, C., and Schoemaker, J. (1998). Increased levels and pulsatility of folliclestimulating hormone in mothers of hereditary dizygotic twins. J Clin Endocrinol Metab 83, 481-486. 585 586 587 51. Van der Stroom, E.M., Konig, T.E., Vink, J.M., Boomsma, D.I., and Lambalk, C.B. (2013). Ovarian reserve and anti-Mullerian hormone (AMH) in mothers of dizygotic twins. Twin Res Hum Genet 16, 634-638. 588 589 52. Hefler, L.A., and Gregg, A.R. (2001). Influence of the angiotensinogen gene on the ovulatory capacity of mice. Fertil Steril 75, 1206-1211. 590 591 592 53. Wide, L., and Eriksson, K. (2013). Dynamic changes in glycosylation and glycan composition of serum FSH and LH during natural ovarian stimulation. Ups J Med Sci 118, 153-164. 25 593 594 595 596 54. Ruth, K.S., Campbell, P.J., Chew, S., Lim, E.M., Hadlow, N., Stuckey, B.G., Brown, S.J., Feenstra, B., Joseph, J., Surdulescu, G.L., et al. (2015). Genome-wide association study with 1000 genomes imputation identifies signals for nine sex hormone-related phenotypes. Eur J Hum Genet. 597 598 599 55. Magill, J.C., Ciccone, N.A., and Kaiser, U.B. (2013). A mathematical model of pulsecoded hormone signal responses in pituitary gonadotroph cells. Math Biosci 246, 3846. 600 601 56. Hillier, S.G. (2009). Paracrine support of ovarian stimulation. Mol Hum Reprod 15, 843850. 602 603 604 57. Laisk-Podar, T., Kaart, T., Peters, M., and Salumets, A. (2015). Genetic variants associated with female reproductive ageing - potential markers for assessing ovarian function and ovarian stimulation outcome. Reprod Biomed Online. 605 606 58. Kirkpatrick, B.W., and Morris, C.A. (2015). A Major Gene for Bovine Ovulation Rate. PLoS One 10, e0129025. 607 608 609 59. Bernard, D.J. (2004). Both SMAD2 and SMAD3 mediate activin-stimulated expression of the follicle-stimulating hormone beta subunit in mouse gonadotrope cells. Mol Endocrinol 18, 606-623. 610 611 60. Gong, X., and McGee, E.A. (2009). Smad3 is required for normal follicular folliclestimulating hormone responsiveness in the mouse. Biol Reprod 81, 730-738. 612 613 614 61. Liu, Y., Chen, X., Xue, X., Shen, C., Shi, C., Dong, J., Zhang, H., Liang, R., Li, S., and Xu, J. (2014). Effects of Smad3 on the proliferation and steroidogenesis in human ovarian luteinized granulosa cells. IUBMB Life 66, 424-437. 615 616 62. Lambalk, C.B., and Boomsma, D.I. (1998). Twinning, cancer, and genetics. Lancet 351, 909-910. 617 618 619 63. Swerdlow, A.J., De Stavola, B.L., Swanwick, M.A., and Maconochie, N.E. (1997). Risks of breast and testicular cancers in young adult twins in England and Wales: evidence on prenatal and genetic aetiology. Lancet 350, 1723-1728. 620 621 622 64. Looyenga, B.D., and Hammer, G.D. (2007). Genetic removal of Smad3 from inhibin-null mice attenuates tumor progression by uncoupling extracellular mitogenic signals from the cell cycle machinery. Mol Endocrinol 21, 2440-2457. 623 Figures titles and legends Figure 1. Manhattan Plot for the genome-wide association meta-analysis of mothers of dizygotic twins. SNPs are plotted on the x-axis according to their position on each chromosome against the significance of the association on the y-axis (shown as –log10 p 26 values). The dotted red line denotes p=5×10–8 statistical significance. The arrows point to the chromosomal region that reached genome-wide significance level. Figure 2. Forest plots depicting risk allele odds ratio (OR) estimates at (A) rs11031006 (near FSHB), (B) rs17293443 (SMAD3) and (C) rs12064669 from each study, the meta-analysis, deCODE (replication) and the overall combined results. Black squares indicate the OR and horizontal lines represent the 95% CIs. The combined results are indicated by the black diamond. Tables titles and legends Table 1. Number of cases and controls included in the meta-analysis and replication; mean maternal age at delivery Table 2. Sample characteristics for fertility measures traits Table 13. Genome-wide significant loci in the meta-analysis of mothers of spontaneous DZ twins (n = 1,980) and replication (n = 3,597) versus screened controls Table 24. Association results for the implicated DZ twinning SNPs in fertility-related measures 27 Table 1. Number of cases and controls included in the meta-analysis and replication; mean maternal age at delivery No. of subjects Cases Controls Mean maternal Cohort age at delivery (SD) Netherlands 806 4,535 29.7 (4.1) Australia 606 6,556 29.8 (4.8) Minnesota 568 1,862 28.3 (4.5) 1,980 12,953 29.3 (4.5) 3,597 (1,356 genotyped; 2,241 297,348 (76,342 genotyped; 30.3 (5.8) in silico) 221,006 in silico) Total discovery Replication Iceland Table 2. Sample characteristics for measure of fertility Traits Sample size Mean (SD) Age at menarche 259,247 -1 Age at menopause 69,360 -1 FSH levels 17,997 - Has children 41,946 0.942 Number of children 41,946 2.9 (1.6) Age at first child 41,946 23.0 (4.7) Age at last child 41,946 32.1 (5.5) Average birth interval 41,946 4.7 (2.5) 5,184 cases / 82,759 controls -1 Polycystic ovary syndrome 1 2 this is a large meta-analysis of summary statistics and sample mean is not available fraction with children; 28 Table 13. Genome-wide significant loci in the meta-analysis of mothers of spontaneous DZ twins (n=1,980) and replication (n=3,597) versus screened controls Meta-analysis SNP Locus Positiona Gene Annotation Risk Replication RAF OR (95%CI) p RAF 0.85 1.41 1.54×10-9 0.85 1.57×10-8 1.23×10-8 OR (95% CI) Combined p p 1.14 (1.06-1.22) 3×10-3 1.25×10-10 0.21 1.15 (1.07-1.23) 1.44×10-4 6.29×10-11 0.07 1.01 (0.89-1.13) 0.88 2.09×10-7 allele rs11031006 11p14.1 30226528 FSHB 5’ upstream G (1.29-1.53) rs17293443 15q22.33 67437863 SMAD3 Intron C 0.24 1.27 (1.19-1.35) rs12064669 1q42.13 230688643 Intergenic C 0.10 1.43 (1.31-1.55) a SNP position according to NCBI Human Genome Build 37; RAF, risk allele frequency; OR, odd ratio; 95% CI, 95% confidence interval 29 Table 24. Association results for the implicated DZ twinning SNPs in fertility-related measures rs11031006-G allele (near FSHB) DZ twinning and fertility-related measure Effect OR (95% CI) / p value rs17293443-C allele (in SMAD3) Effect Beta (95% CI) OR (95% CI) / p value Beta (95% CI) DZ twinninga Increase 1.41 (1.29, 1.53) 1.54×10-9 Increase 1.27 (1.19, 1.35) 1.57×10-8 FSH levels (SD units)b Increase 0.11 (0.078, 2.3×10-10 ̶ 0.016 (-0.014, 0.30 0.15) Age at menarche (years)c,d Earlier -0.04 (-0.012, 0.045) 8.5×10-10 Earlier̶ 0.011) Age at menopause (years)e Earlier -0.2165 (-0.065, ̶ 8.5×10-14 0.009 (-0.048, 0.049) 1.07 (0.98, 1.17) 0.12 ̶ 0.96 (0.89, 1.04) 0.34 Increase 0.014 (0.00091, 0.03 ̶ 0.0048 (-0.0062, 0.39 0.27)h Earlier -0.20 (-0.31, - 0.016)h ̶ 5.3×10-4 0.032 (-0.065, 0.086) Age at last child (years)f Earlier -0.097 (-0.21, ̶ 0.015 (-0.037, 0.51 0.13) 0.08 Later 0.015) Average birth interval (years)f ‒0.71 ̶ Has children (yes/no)f Age at first child (years)f 0.84 0.010) 0.052) Number of childrenf -0.001 (-0.127, 0.14 (0.043, 4.7×10-3 0.24) 0.57 ̶ 0.0051 (-0.040, 0.82 30 0.067) PCOSg Decrease 0.90 (0.84, 0.95) 0.050) 12.39×10- ̶ 1.00 (0.96,1.05) 0.78 49 MODZT GWAS; bdeCODE sample of 17,997 individuals with FSH levels; c,dLunetta dDay et all . and Perry et al. (ref 21 and 2221) and Perry et al. (ref 22); eDay et al. (ref 2424) and Stolk et al. (ref 25); fdeCODE sample of 41,946 women; gPolycystic ovary syndrome, Hayes Day et al. (ref 223); h factor of increase from log-linear model.. a 31