{"paper_id":"59ba4052-94e0-4899-9b01-e48dbf6d6f52","body_text":"Human reproductive behavior – age at first birth (AFB) and number of children ever born (NEB) – has been associated with human development, 1 , 2  infertility 3 , 4  and neuropsychiatric disorders 5 . Reproductive tempo (AFB) and quantum (NEB) are cross-cutting topics in the medical, biological, evolutionary and social sciences and central in national and international policies. 6  Advanced societies experienced a rapid postponement of AFB, with the mean AFB of women now being 28–29 years in many countries. 7  This increase in AFB has been linked to lower fertility rates, unprecedented childlessness (~20%) and infertility, which affects 10 to 15 % of couples. 8  An estimated 48.5 million couples worldwide are infertile, with a large part of subfertility, particularly in men, remaining unexplained. 9  Although infertility has been related to advanced AFB, ovulation defects, spermatogenic failure, and single- or polygenic defects, their causal effects remain unsubstantiated. 10\nUntil now, genetic and clinical research has focussed on proximal infertility phenotypes. 3 , 4 , 10 , 11  AFB and NEB represent accurate measures of complex reproductive outcomes, are frequently recorded and consistently measured, and are key parameters for demographic population forecasting. 12  There is evidence of a genetic component underlying reproduction, with heritability estimates of up to 50% for AFB and NEB ( Supplementary Figure 1 ). 6  A recent study attributed 15% of the variance of AFB and 10% of NEB to common genetic variants. 13  There are also sex-specific differences in human reproduction, related to the timing of fertility, fecundability and sex-genotype interactions ( Supplementary Note ). Researchers have given less attention to traits such as NEB due to an erroneous and frequently repeated misinterpretation of Fisher’s 14  Fundamental Theorem of Natural Selection that the additive genetic variance in fitness should be close to zero. The misreading had a naively intuitive appeal: genes that reduce fitness should have been less frequently passed on. Fisher, however, actually argues that fitness is moderately heritable in human populations (for a discussion see the  Supplementary Note ). Since no established genes are currently available for clinical testing of infertility, 10  isolating genetic loci and their causative effects has the potential to provide new insights into the etiology of reproductive diseases and novel diagnostic and clinical technologies for infertility treatment.\n\nWe report the largest meta-analysis of genome-wide association studies (GWAS) to date of 251,151 individuals for AFB and 343,072 for NEB from a total of 62 cohorts of European ancestry. We identify 12 independent loci (10 of which are novel and 2 previously identified in a study on age at first sexual intercourse 11 ) that are significantly associated with AFB and/or NEB in men, women and/or both sexes combined ( Table 1 ). Follow-up analyses identified a number of genetic variants and genes that likely drive GWAS associations. We also quantified the genetic overlap with biologically adjacent reproductive, developmental, neuropsychiatric and behavioral phenotypes. A detailed description of all materials and methods is available in the  Supplementary Note .\nAssociations of AFB (mean ± SD 26.8±4.78 years) and/or NEB (mean ± SD 2.3±1.43 children) with NCBI build 37 HapMap Phase 2 imputed SNPs were examined in 62 cohorts using multiple linear regression under an additive model, in men and women separately ( Supplementary Note ). Associations were adjusted for principal components – to reduce confounding by population stratification 15  – as well as for the birth year of the respondent and its square and cubic to control for non-linear birth cohort effects ( Supplementary Note  and  Supplementary Tables 1–,2  ). NEB was assessed only for those who had completed their reproductive period (age ≥ 45 women; age ≥ 55 men), while AFB was only assessed for those who were parous. Quality control (QC) was conducted in two independent centers using QCGWAS 16  and EasyQC 17  ( Supplementary Note ). Results were subsequently meta-analyzed for the 2.4M SNPs that passed QC filters ( Supplementary Note ) and reported for men and women combined and separately.\nWe applied a single genomic control at the cohort level and meta-analyzed results using a sample-size weighted fixed effect method in METAL ( Supplementary Note ). The PLINK clumping function isolated ‘lead SNPs’ – i.e. those with the lowest  P -value for association that are independently associated – using an r 2  threshold of 0.1 and distance threshold of 500 kb. For AFB, we identified ten genome-wide significantly associated loci (i.e.,  P <5x10 −8  for combined and  P <1.67x10 −8  for sex-specific results adjusted for multiple testing) of which nine were significantly associated in both sexes combined and one in women only (N=154,839) ( Figure 1a ,  Table 1 ). Three loci were significantly associated with NEB: two in both sexes combined and one in men only (N=103,736) ( Figure 1b ,  Table 1 ,  Supplementary Note ). One locus on Chr 1 reached significance for association with both AFB and NEB with an r 2  of 0.57 between the two lead SNPs, suggesting a shared genetic basis for the two traits ( Table 2 ). A statistical test of sex-specific effects confirms that differences are mainly due to variation in sample size and not variation in effect sizes ( Supplementary Note ).\nAs for other complex traits 18 , the Q-Q plots of the meta-analyses exhibit strong inflation of low  P -values ( Figure 2 ), suggesting that although cohorts controlled for the top principal components and cohort-level genomic control was applied ( Supplementary Note ), residual population stratification may remain. However, the LD Score intercept method 19  as well as a series of individual and within-family regression analyses using polygenic scores as predictors 20 , 21  ( Supplementary Note ) indicated that the observed inflation is almost entirely attributable to a true polygenic signal, rather than population stratification.\nTo increase the power to find statistically significant associations and causal genes, we additionally performed a gene-based GWAS using VEGAS. 22 , 23  The results confirmed top hits from the single-SNP analyses. For AFB, seven loci from the SNP-based GWAS were also represented in the gene-based analysis ( Supplementary Table 3 ), and three additional loci emerged, represented by  SLF2  (Chr 10),  ENO4  (Chr 10) and  TRAF3-AMN  (Chr 14). For NEB, one locus from the SNP-based GWAS was represented in the gene-based analysis – i.e.  GATAD2B  (Chr 1) – and one novel locus on Chr 17 was identified ( Supplementary Table 4 ).\nTo identify functional and potentially causal variants – coding or regulatory – within loci identified in the SNP-based GWAS ( Table 1 ), we first performed an  in silico  sequencing annotation analysis using the post-GWAS pipeline reported by Vaez  et al . 24  This showed that rs10908557 on Chr 1 is in high LD with non-synonymous SNPs in  CRTC2  (rs11264680; r 2 =0.98) and  CREB3L4  (rs11264743; r 2 =0.94) (see Causal genes,  Supplementary Table 5 ). Interestingly, rs11264743 is considered ‘deleterious’ and ‘probably damaging’ by SIFT and PolyPhen, respectively (Ensembl release 83). In addition, rs2777888 on Chr 3 is in high LD with two non-synonymous SNPs in  MST1R  (rs2230590; r 2 =0.95 and rs1062633; r 2 =0.95) ( Table 1 ,  Supplementary Table 5 ).\nWe subsequently performed a comprehensive analysis using results from the ENCODE 25  and RoadMap Epigenomics 26  projects as integrated in RegulomeDB, 27  to identify variants that likely influence downstream gene expression via regulatory pathways. Amongst all SNPs that reached  P <5x10 −8  in the meta-analyses (N=322), 50 SNPs in five loci show the most evidence of having functional consequences ( Table 1 ,  Supplementary Table 6 ). Two sets of SNPs on Chr 1 (18 SNPs) and Chr 3 (25 SNPs) particularly stand out. The most promising SNP in the Chr 1 locus (rs6680140) is located in an H3K27ac mark, often found near active regulatory elements, and lies in a DNaseI hypersensitivity cluster where eight proteins are anticipated to bind. One of these proteins is cAMP responsive element binding (CREB) binding protein, encoded by  CREBBP  (see Causal genes). In the Chr 3 locus, rs2526397 is located in a transcription factor-binding site and is an eQTL for  HYAL3  in monocytes, while rs2247510 and rs1800688 are located in H3K27ac sites and DNaseI hypersensitivity clusters where a large number of transcription factors are expected to bind (see Causal genes,  Supplementary Table 6 ). An analysis using Haplotter showed that rs2247510 and rs7628058 in the Chr 3 locus are amongst the 5% of signals that show most evidence of positive selection in the population. The same applies to the lead SNP of the Chr 14 locus for NEB (rs2415984).\nInformation on the function and anticipated relevance of genes in the 12 loci identified in the SNP-based GWAS that are most likely to be causal based on all evidence discussed below is provided in  Table 2 .\nIdentifying alterations in gene methylation status and/or expression levels in relation to GWAS-identified variants may help prioritize causal genes. We examined associations with gene expression and methylation status for the 12 independent lead SNPs in whole-blood BIOS expression (e)QTL (N=2,116) and methylation (me)QTL databases in  cis  and  trans  (N=3,841).  28 , 29  Seven SNPs were associated with gene expression in  cis  of 54 unique genes ( Table 1 ,  Supplementary Table 7 ). Five of the seven SNPs were in high LD ( r 2 >0.8) with the strongest eQTL of at least one of the genes within the corresponding loci, indicating that the SNP associated with AFB or NEB and the SNP most significantly associated with expression tag the same functional site, i.e., rs10908557 (associated with the expression of  CRTC2  and  SLC39A1 ), rs1160544 ( AFF3 ), rs2777888 ( RBM6 ,  RNF123  and  RBM5 ), rs2721195 ( CYHR1 ,  GPT ,  RECQL4  and  PPP1R16A ) and rs293566 ( NOL4L ). Three SNPs were associated with the expression of a total of eight genes in  trans  ( Table 1 ,  Supplementary Table 8 ). Of these SNPs, only rs2777888 was in high LD ( r 2 >0.8) with the strongest eQTL for three of its five associated genes:  LRFN1 ,  LAMP2  and  FGD3 .\nThe meQTL analysis showed that 11 of the 12 independent lead SNPs were associated with DNA methylation of a total of 131 unique genes in  cis  ( Table 1 ,  Supplementary Table 9 ). Seven of the 11 SNPs were in high LD ( r 2 >0.8) with the strongest meQTL of one of the corresponding methylation sites, i.e., rs10908557 (associated with methylation of  CRTC2, SLC39A1, CREB3L4, DENND4B  and  RAB13 ), rs1160544 ( AFF3 ), rs2777888 ( CAMKV ), rs6885307 ( C5orf34 ), rs10056247 ( JADE2 ), rs2721195 ( CYHR1 ) and rs13161115 ( EFNA5 ). Three of the SNPs were associated with the same genes for both methylation and gene expression in  cis : rs10908557 ( CRTC2 ), rs1160544 ( AFF3 ) and rs2721195 ( CYHR1 ) ( Supplementary Tables 7,9 ). Three SNPs were associated with methylation of a total of ten genes in  trans  ( Table 1 ,  Supplementary Table 10 ). Of these SNPs, only rs2777888 was in high LD ( r 2 >0.8) with the strongest meQTL of a corresponding methylation site ( ASAP3 ). Of note: rs2777888 was also a  trans  eQTL.\nWe used four publicly available bioinformatics tools to systematically identify genes that are more likely than neighboring genes to cause the associations identified by our GWAS. Of all genes that reached  P <0.05 in analyses using Endeavour, 30  MetaRanker 31  and ToppGene, 32  eight genes were prioritized for both AFB and NEB:  TPM3 ,  GRM7 ,  TKT ,  MAGI2 ,  PTPRD ,  PTPRM ,  RORA  and  WT1 . Data-driven Expression-Prioritized Integration for Complex Traits (DEPICT) – a fourth, comprehensive and unbiased recently described gene prioritization tool 33  – identified three genes in GWAS significant loci as likely being causal for AFB ( MON1A ,  RBM6  and  U73166.2 ) ( Supplementary Tables 11, 12 ).\nAn analysis using RegulomeDB identified 50 SNPs in five loci that most likely have regulatory consequences (see Causal variants,  Supplementary Table 6 ). Eighteen and 25 of these SNPs are within the Chr 1 and Chr 3 loci, respectively. Amongst the genes that – at a protein level – bind at the site of one or more of the 18 variants in the locus on Chr 1 are  CREBBP ,  HNF4A ,  CDX2  and  ERG . These genes may act upstream in the causal pathway and influence the expression of causal genes at this locus. Of the 25 SNPs on Chr 3, ten were eQTLs for  RBM6  in monocytes, and seven were eQTLs for  HYAL3  in monocytes. Amongst the genes that – at a protein level – bind at rs2247510 and rs1800688 in the Chr 3 locus are  ARID3A ,  REST  and  TFAP2C , as well as  HNF4A  and  CDX2 , which also bind at the Chr 1 locus.\nFive genes encode proteins that bind at the site of both SNPs on Chr 2 that reach  P <5x10 −8  in the meta-analysis of GWAS. One of these is  REST  (see Chr 3 locus), another one –  ESR1  – is the most likely causal gene in the Chr 6 locus.\nFunctional network analysis using five prioritized candidate gene sets as input ( Supplementary Note ) showed no significantly enriched biological function. No biological function was significantly enriched after correction for multiple testing using the Benjamini-Hochberg procedure. Similarly, no reconstituted gene sets and cell or tissue types were significantly enriched in the GWAS meta-analysis results based on results from the DEPICT analysis ( Supplementary Tables 13–20 ). The lack of significantly enriched genes, tissue sets and biological functions reflects the need for a larger sample size but also the distal nature of the phenotypes, which are influenced by a mixture of biological, psychological and socio-environmental factors.\nTo assess the predictive power of our results, we produced polygenic scores for AFB and NEB using sets of SNPs whose nominal  P -values ranged from  P <5x10 −8  (i.e. using only genome-wide significant SNPs) to 1 (using all SNPs that passed quality control) using PRSice 34  ( Supplementary Note ). We then performed a series of four different out-of-sample predictions in four independent cohorts: HRS, Lifelines, STR and TwinsUK. Across the four cohorts, the mean predictive power of a polygenic score constructed from all measured SNPs is 0.9% for AFB and 0.2% for NEB ( Supplementary Figure 2 ). Despite the low predictive power of the polygenic scores, the results showed that a 1 standard deviation (SD) increase of the NEB polygenic score is associated with a 9% (95% CI 5%–13%) decrease in the probability for women to remain childless, with no significant association in men ( Supplementary Table 21 ). When we control for right-censored data using a survival model for AFB, we found that a 1SD increase in the AFB polygenic score is associated with an 8% (95% CI 7%–10%) reduction in the hazard ratio of reproduction in women and 3% (95% CI 1%–5%) in men ( Supplementary Table 22 ). As an additional test, we examined whether the aforementioned polygenic scores for AFB and NEB can predict related fertility traits such as age at menopause and age at menarche ( Supplementary Table 23 ). Our estimates indicate that a 1SD increase of the AFB polygenic score is associated with a 3% decrease in age at natural menopause (95% CI 1%–5% ) and a 20 day increase in age at menarche (95% CI 0.4–40 days).\nSeveral loci for which the associations with AFB and NEB reach genome-wide significance are associated with behavioral and reproductive phenotypes. The lead SNPs in the Chr 2 and Chr 3 loci have been associated with educational attainment 35  and the locus on Chr 5 with age at menarche 36  while the locus on Chr 6 has recently been associated with age at first sexual intercourse 37  ( Supplementary Table 24 ). Some of the 38 loci for age at first sexual intercourse that were recently identified in 125,667 UK Biobank participants were also associated with AFB (in/near  RBM6 – SEMA3F  and  ESR1 ) and NEB (in/near  CADM2  and  ESR1 ). The lead SNPs for  RBM6 – SEMA3F  (rs2188151) and  ESR1  (rs67229052) are in LD with our lead SNPs for AFB on Chr 3 (r 2 = 0.44) and Chr 6 (r 2 =0.94), respectively. An  in silico  pleiotropy analysis showed that our lead SNP in the Chr 3 locus (rs2777888) is in LD (r 2 =0.59) with rs6762477 – which has been associated with age at menarche 2  – while other SNPs in the same locus have been associated with HDL cholesterol 38  (rs2013208; r 2 =0.81) and BMI 39  (rs7613875; r 2 =0.81) ( Supplementary Table 5 ). We next performed an exploratory analysis using the proxy-phenotype method 40  to examine if the SNPs strongly associated with AFB in women are empirically plausible candidate SNPs for related traits like age at menarche and age at menopause ( Supplementary Note ). After controlling for multiple testing, we identified three AFB-associated SNPs near rs2777888 on Chr 3 (rs9589, rs6803222 and rs9858889) that are also associated with age at menarche ( P <4.10x10 −4 ). None of the AFB or NEB-associated SNPs are associated with age at menopause.\nWe performed a bivariate LD score regression analysis 41  to estimate the pairwise genetic correlation with 27 publicly available GWAS results for traits associated with human reproductive behavior ( Supplementary Note ). AFB shows significant and positive genetic correlations with the human (reproductive) developmental traits age at menarche, voice breaking, age at menopause, birth weight and age at first sexual intercourse, as well as with years of education. Conversely, having more AFB-increasing alleles is associated with a lower genetic risk of smoking (ever, number of cigarettes, later onset) and with lower insulin resistance-related phenotypes, i.e. BMI, waist-hip-ratio adjusted for BMI, fasting insulin, triglyceride levels and risk of diabetes ( Figure 3  and  Supplementary Table 25 ). All genetic correlations remain significant after Bonferroni correction for multiple testing ( P <2.6x10 −3 ). Years of education ( P =6.6x10 −14 ) and age at first sexual intercourse ( P =1.14x10 −15 ) are the only traits that show significant and negative genetic correlations with NEB. We also observed significant genetic correlations of  r g =0.86 (SE=0.052) for AFB and  r g =0.97 (SE=0.095) for NEB between men and women, implying that most genetic effects on reproductive behavior resulting from common SNPs are shared across the sexes.\n\nThis GWAS is the largest genetic epidemiological discovery effort for human reproduction to date, with critical implications for population fitness and clear physiological mechanisms linking hypothesized genes and observed phenotypes. Related studies previously focussed on reproductive life span 42 , 43 , age at first sexual intercourse 11  and more proximal infertility phenotypes, 2 – 4  largely overlooking AFB and NEB. The rapid postponement of AFB and increased infertility and involuntary childlessness in many societies 7  makes it important to uncover the genetic and biological architecture of reproduction. We identify ten novel and confirm two recently identified genetic loci that are robustly associated with AFB and NEB, as well as variants and genes within these loci that likely drive these associations. Four additional loci were identified in a gene-based GWAS.\nTwo loci that show interesting results in follow-up analyses are located on Chrs 1 and 3. The lead SNPs of the Chr 1 locus for AFB and NEB are in LD with likely functional non-synonymous SNPs in genes encoding: 1) CREB (cAMP responsive element binding) regulated transcription co-activator 2 ( CRTC2 ), which at a protein level acts as a critical signal mediator in follicle-stimulating hormone (FSH) and transforming growth factor β1(TGFβ1)-stimulated steroidogenesis in ovarian granulosa cells 44 ; and 2) CREB protein 3-like 4 ( CREB3L4 ), 45  which in humans is highly expressed in the prostate, ovaries, uterus, placenta and testis, and plays a role in spermatid differentiation 46  and male germ cell development. 47  The lead SNP and/or functional variants in LD with it are also associated with the methylation status of these two genes and expression of  CRTC2  in whole blood and lymphocytes. Three promising functional variants in the Chr 1 locus reside in binding sites of the transcriptional co-activator CREB binding protein (CREBBP). In addition to a direct effect of the above-mentioned non-synonymous SNPs on protein function, the associations of AFB and NEB with variants in the locus on Chr 1 may thus be mediated by alterations in cAMP responsive element binding in men and women. The locus on Chr 1 also harbours  DENND4B , a paralogue of  DENND1A , implicated in PCOS. 48  While  DENND1A  is expressed at the protein level in the ovary and testis, DENND4B is in the cervix, and its function and role are less well understood.\nThe lead SNP of the locus on Chr 3 (rs2777888) is associated with methylation and expression of several genes – in  cis  and  trans  – that are known to play a role in cell cycle progression and/or sperm function. First, rs2777888 is associated with the expression of  RNF123  (or  KPC1 ) in  cis , which plays a role in cellular transition from the quiescence to proliferative state. Secondly, rs2777888 – or functional variants in LD with it – may influence the cell cycle by altering the expression of  RBM5  and  RBM6  (RNA binding motif proteins 5 and 6). The former plays a role in cell cycle arrest and apoptosis induction and regulates haploid male germ cell pre-mRNA splicing and fertility in mice.  RMB5  mutant mice exhibit spermatid differentiation arrest, germ cell sloughing and apoptosis, leading to lack of sperm in ejaculation. 49  Thirdly, rs2777888 affects expression of  LAMP2  in  trans , which is located on the X chromosome and encodes a lysosomal membrane protein involved in the acrosome reaction, i.e. the enzymatic drill allowing sperm to penetrate and fertilize ovum. 50 \n LAMP2  is expressed at the protein level in male and female reproductive organs with an effect size of rs2777888 for  LAMP2  mRNA expression that is almost twice as large in women than in men ( Supplementary Figure 4 ). This suggests an important role for the lysosome in both sperm and ovum. Finally, functional variants in the Chr 3 locus are associated with the mRNA expression of  HYAL3  in monocytes (hyaluronoglucosaminidase 3). The latter degrades hyaluronan, which also plays an important role in sperm function and the acrosome reaction. 49 , 51\nFunctional follow-up experiments could disentangle the potential interplay between many candidate genes in the loci on Chrs 1 and 3 on reproductive behavior and – given our  in silico  results – infertility. This can be extended to candidate genes in the remaining loci identified in the present study, some of which are relevant for fertility: mice lacking  EFNA5  (Chr 5 NEB locus) are subfertile, 52 \n ESR1  on Chr 6 encodes an estrogen receptor,  53 , 54  and  CYHR1  on Chr 8 is involved in spermatogenesis 55 . Such experiments would help understand whether binding of estrogen receptor 1 – encoded by  ESR1  in the locus on Chr 6 – at the site of functional variants in the locus on Chr 2 drives or mediates the association with AFB in the Chr 2 locus, as well as to identify and characterize causal genes. Recent developments in high-throughput, multiplex mutagenesis using Clustered Regularly Interspaced Short Palindromic Repeats (CRISPR) and associated systems (Cas9) allow such experiments to be performed using  in vivo  model systems. 56\nAFB and NEB are not only driven by biological processes, but are also subject to individual choice and personal characteristics – such as personality traits – as well as by the historical, cultural, economic and social environment (e.g., effective contraception, childcare availability). Demographic research has shown a strong positive association between AFB and educational attainment. 12  We show that the associations between fecundity, reproductive behavior and educational attainment are partly driven by genetic factors, and identified loci that are associated with AFB as well as with e.g., age at first sexual intercourse  37  and educational attainment. 35\nOur findings are anticipated to lead to insights into how postponing reproduction may be more detrimental for some – based on their genetic make-up – than others, fuel experiments to determine ‘how late can you wait’ 57  and stimulate reproductive awareness. Causal genes in the loci we identified may serve as novel drug targets, to prevent or delay age-related declines in fertility and sperm quality, or increase Assisted Reproductive Technology efficiency. Our study is the first to examine the genetics of reproductive behavior in both men and women, and the first that is adequately powered to identify loci both in women and men. We also provide support for Fisher’s theorem that fitness is moderately heritable in human populations. While effect sizes of the identified common variants are small, there are examples of GWAS-identified loci of a small effect that end up leading to important biological insights. 58 , 59  Many of our findings suggest a role for sperm quality, which is one lead for researchers to pursue. Since current sperm tests remain rudimentary, our findings could serve as a basis for ‘good quality’ sperm markers. We identified variants that are likely causal – both coding and regulatory – as well as a set of genes that likely underlie the associations we identified. Follow-up experiments in animal models are required to confirm and characterize the causal genes in these loci.\n\nGenome-wide association analyses of age at first birth (AFB) and number of children ever born (NEB) were performed at the cohort level according to a pre-specified analysis plan (see  Supplementary Note ). Cohorts uploaded results imputed using the HapMap 2 CEU (r22.b36) or 1000G reference sample. Cohorts were asked to only include participants of European ancestry, with no missing values on all relevant covariates (sex, birth year, and cohort specific covariates), who were successfully genotyped genome-wide, and passed cohort-specific quality controls. We followed the QC protocol of the GIANT consortium’s recent study of human height 87  and employed QCGWAS 88  and EasyQC 89  software, which allowed us to harmonize the files and identify possible sources of errors in association results.\nCohort association results (after applying the QC filters) were combined using sample-size weighted meta-analysis with genomic control (GC) correction within each study, implemented in METAL. 90  SNPs were considered genome-wide significant at  P -values smaller than 5×10 −08  (α of 5%, Bonferroni-corrected for a million tests). The meta-analyses were carried out by two independent analysts. Detailed results of each genome-wide significant locus are shown in in  Supplementary Figures 4–29 .\nThe total sample size of the meta-analysis is  N =251,151 for AFB pooled and  N =343,072 for NEB pooled. The PLINK clumping function 91  was used to identify the most significant SNPs in associated regions (termed “lead SNPs”). Detailed cohort descriptions, information about cohort-level genotyping and imputation procedures, cohort-level measures, and quality control filters are shown in  Supplementary Tables 26, 27  and discussed in the  Supplementary Note .\nWe applied a method recently developed by Zhu and colleagues 92  to estimate dominant genetic effects based on the genetic relatedness of unrelated individuals. Our results based on the combined samples of TwinsUK and Lifelines show no evidence for dominant genetic effects for either NEB (1.0x10 −07 , SE=0.07,  P =0.45) nor AFB (0.02, SE=0.08,  P =0.43. Results are shown in  Supplementary Table 28  and discussed in the  Supplementary Note .\nAs joint analysis of correlated traits may boost power for mapping functional loci, we applied a recently developed multi-trait analysis method 93  to test the association between each variant and the two correlated traits AFB and NEB simultaneously using multivariate analysis of variance (MANOVA) (see  Supplementary Note  and  Supplementary Table 29 ). The analysis was performed based on the genome-wide meta-analysis summary statistics of each single trait. Although it did not reveal additional genome-wide significant loci (λ=0.995), it accounted for the correlation between the two phenotypes, thus improving the strength of two signals on chromosomes 1 and 5, indicating possible pleiotropic architecture between AFB and NEB ( Supplementary Figure 30 ). The analysis also provided a conditional association test of the genetic effect of each variant on AFB including NEB as a covariate, and on NEB including AFB as a covariate ( Supplementary Figure 31 )\nWe used two methods to assess whether our GWAS results exhibited signs of population stratification (see  Supplementary Note ). First, we used the LD Score intercept method described in Bulik-Sullivan  et al . 94  to test whether inflation in chi-square statistics was due to confounding biases such as cryptic relatedness and population stratification. In all six cases, the intercept estimates were not significantly different from 1, suggesting no appreciable inflation of the test statistics attributable to population stratification. Second, we conducted a series of individual and within-family (WF) regressions using polygenic scores (PGS) as predictors 95 – 97  on a dataset of DZ twins (STR and TwinsUK). The regression analyses showed that WF regression coefficients for both AFB and NEB were statistically different from zero when the  P -value threshold was sufficiently high ( Supplementary Tables 30, 31  and  Supplementary Figures 32, 33 ).\nIn addition to the pooled GWAS results presented in the main text, we also ran sex-specific GWAS meta-analyses for AFB and NEB (see  Supplementary Note ). The sample size for sex-specific analysis was: AFB women,  N =189,656; AFB men,  N =48,408; NEB women  N =225,230; NEB men  N =103,909. Our results indicated 6 genome-wide significant ( P <5x10 − 08 ) independent SNPs for AFB women and 1 genome-wide significant independent SNP for NEB men ( Supplementary Table 32  and  Supplementary Figures 34, 35 ). We also used LD score bivariate regression and GREML bivariate analysis to estimate the genetic correlation among men and women based on the sex-specific summary statistics of AFB and NEB meta-analysis results. Our estimates based on LD bivariate regression indicated a genetic correlation of r g =0.86 (SE=0.052) among the sexes for AFB and r g =0.97 (SE=0.095) for NEB. Results are shown in  Supplementary Tables 33, 34  and discussed in the  Supplementary Note .\nWe performed out-of-sample prediction and calculated polygenic scores for AFB and NEB, based on GWA meta-analysis results and used regression models to predict the same phenotypes in four independent cohorts: HRS, Lifelines, STR and TwinsUK (see  Supplementary Note  and  Supplementary Figure 2 ). We ran ordinary least-squares (OLS) regression models and reported the R 2  as a measure of goodness-of-fit of the model. In addition, we tested how well our polygenic scores for NEB could predict childlessness at the end of the reproductive period (using age 45 for women and 55 for men),  Supplementary Table 21 . Since age at first birth is observed only in parous women, we adopted an additional statistical model to account for censoring (Cox Proportional hazard model,  Supplementary Table 22 ) and selection (Heckman selection model,  Supplementary Table 35 ). We additionally tested the predictive value of our polygenic scores for AFB for age at menarche (using TwinsUK) and age at menopause (using Lifelines),  Supplementary Table 23 . Finally, we examined whether menopause variants are associated with AFB. We calculated a PGS of age at menopause based on the recent GWAS results from Day et al. (2015) 98  and applied them to LifeLines and TwinsUK ( Supplementary Table 36 ).\nWe used information from 27 publicly-available GWAS results to estimate the amount of genetic correlations between AFB and NEB and related traits ( Supplementary Table 25  and  Figure 3  in the main text) using LD score bivariate regression (see  Supplementary Note ). Details on these phenotypes are provided in the  Supplementary Note . A conservative Bonferroni-corrected  P -value threshold of  P <1.85×10 −03  (=0.05/27) was used to define significant associations. We also tested the correlation between NEB and AFB using a bivariate GREML analysis on the Women’s General Health Study (WGHS,  N =40,621).\nFollowing the results on genetic overlap with other phenotypes we tested – in a quasi-phenotype replication setting – whether the SNPs strongly associated with AFB in women were empirically plausible candidate SNPs for age at menarche and age at menopause (see  Supplementary Note ). We used a two-stage approach applied in other contexts. 99 , 100  In the first stage, we conducted a meta-analysis of AFB excluding the cohorts that were part of the meta-analysis of the phenotype we intended to replicate. We merged the SNPs from this meta-analysis with the publically available association results of the most recent GWAS on age at menarche 2  and age at menopause 101  from the Reprogen consortium website. 1  SNPs that were not present in both files were dropped from the analysis. We aligned the alleles and the direction of effects using the EasyStrata software. 102  We then selected the independent SNPs with a  P -value<1x10 −05 , using the clump procedure in PLINK 91 , (1000Kb window and LD threshold of R 2 >0.1) to identify the most significant SNPs in associated regions included in both files. We defined “prioritized SNP associations” as those that passed the Bonferroni correction for the number of SNPs tested (0.05/122=4.10x10 −04 , both in age at menarche and age at menopause). Our results identified three SNPs after Bonferroni-correction that can be used as good candidates for age at menarche. We did not isolate any clear “candidate SNP” for age at menopause ( Supplementary Figure 36 ).\nWe performed gene-based testing with the full GWAS set (~2.5M HapMap-imputed SNPs) of both phenotypes using VEGAS (see  Supplementary Note  and  Supplementary Tables 3,4 ). 23 , 103  This software has the advantage of accounting for LD structure and the possibility to define a range beyond the gene bounds to include intergenic regions in the analysis. We defined a 50kb extra window surrounding the genes and considered every SNP for the gene-based analysis, ran the analyses per chromosome with up to 10 06  permutations and considered  P <2.5x10 −06  as the threshold for significance (0.05/~20.000 genes).\nFor each of the 12 SNPs identified in the GWAS, local ( cis , exons/methylation sites < 1 MB from the SNP) and genome-wide ( trans , exons/methylation sites > 5 MB from the SNP) effects were identified by computing Spearman rank correlations between SNPs and local or global exons/methylation sites (see  Supplementary Note ). Bonferroni multiple testing correction was performed for the 12 SNPs tested ( P <2.5x10 −06  for  cis  meQTL analysis,  P <1x10 −08  for  trans  meQTL analysis,  P <1.2 x10 −06  for  cis  eQTL analysis,  P <1.3x10 −08  for  trans  eQTL analysis). For each of the significant associations, the exons/methylation sites were selected, the strongest eQTLs were identified for these exons/methylation sites, and LD between the strongest eQTLs and the corresponding SNP identified in the GWAS were computed. LD was computed using BIOS genotypes (the genotypes used for eQTL and meQTL mapping).\nWe used RegulomeDB 27  to identify variants amongst the 322 SNPs that reached  P <5x10 −08  for association with AFB and/or NEB in the meta-analysis of GWAS that likely influenced regulation of gene expression (see  Supplementary Note ). RegulomeDB integrates results from RoadMap Epigenomics 26  and the ENCODE project. 104  SNPs showing the most evidence of being functional – defined as a RegulomeDB score <4 – were subsequently examined in more detail in terms of effects on gene expression (eQTLs) and their protein-binding capacity ( Supplementary Table 6 ).\nPotentially causal genes for the associations identified by GWAS were identified using four previously described bioinformatics tools: ToppGene 4 , Endeavour 5 , MetaRanker 6 , and DEPICT 7 . To this end, we first retrieved positional coordinates for all lead SNPs according to GRCh37/hg19 using Ensembl’s BioMart. These coordinates were used to extract all genes located within ±40kb of lead SNPs using the UCSC table browser. The identified genes then served as input for ToppGene and Endeavour. Genes with established roles in fertility served as training genes in this procedure, i.e.  BRCA1 ,  EGFR ,  ERBB2-4 ,  HSD17B1 ,  RBM5 ,  ESR1 ,  ESR2  and  FSHB . For MetaRanker we provided SNPs that reached  P <5x10 −04  and their chromosomal position as input, together with the previously mentioned set of training genes. Since ToppGene, Endeavour and MetaRanker are biased towards larger and well-described genes, we also performed a gene prioritization procedure using DEPICT. 7  All SNPs that reached  P <5x10 −04  in the meta-analysis served as input, and information on prioritized genes, gene set enrichment, and tissue/cell type enrichment were extracted. Genes were subsequently prioritized that: 1) reached  P <0.05 in DEPICT; or 2) reached  P <0.05 in ToppGene, Endeavour and MetaRanker ( Supplementary Table 37 ).\nDEPICT was used to identify gene set, cell type and tissue enrichment analyses, using the GWAS-identified SNPs with  P <5x10 −04  as input (see  Supplementary Note ). Due to the relatively small number of identified loci, DEPICT was only able to perform these analyses for AFB and NEB pooled, and AFB women. To construct a functional association network, we combined five prioritized candidate gene sets into a single query gene set which was then used as input for the functional network analysis. 24  We applied the GeneMANIA algorithm together with its large set of accompanying functional association data. 105  We used the Cytoscape software platform, 106  extended by the GeneMANIA plugin (Data Version: 8/12/2014, accessed April 24, 2016). 107  All the genes in the composite network, either from the query or the resulting gene sets, were then used for functional enrichment analysis against Gene Ontology terms (GO terms) 108  to identify the most relevant GO terms using the same plugin. 107\nPrevious research based on twin studies shows differential heritability of fertility behavior across birth cohorts. 109 , 110  We used the Swedish Twin Register (STR) to examine whether the effect of a polygenic score (PGS) of AFB and NEB varies across birth cohort. We followed the analysis presented in the recent GWAS of education 111  and divide the sample into six groups based on their year of birth. Each group spans five birth years, with the oldest ranging from 1929–1933 and the youngest born between 1954–1958.  Supplementary Table 38  reports the estimated coefficient from these regressions. The results indicate a U-shaped trend in AFB and a linear decline in NEB, but do not provide any clear evidence of interaction effects between the PGS’s and birth cohort. We additionally tested the interaction effects between educational level and the PGS of AFB and NEB in three different samples (LifeLines, STR and HRS).  Supplementary Table 39  reports the estimated coefficient from these regressions. The results indicate that years of education are positively associated with AFB in both LifeLines and STR, and negatively associated with NEB in the HRS. With the exception of NEB in the HRS, we found no evidence of GxE effects with education.\nTo estimate the robustness of our results for AFB, we conducted two additional analyses. First, we estimated how the coefficients change if we control for Educational Attainment (EA). Using data from deCODE, we ran an additional association analysis using the 10 loci that were genome-wide significant in the meta-analysis ( P <5x10 −08 ). The analysis has been restricted to individuals born between 1910 and 1975, who also had data available on completed education. The total sample size is 42,187 (17,996 males and 24,191 females). The analysis is adjusted for sex, year of birth (linear, squared and cubic), interaction between sex and year of birth and the first 10 PCAs. Education is measured by years of education, ranging between 10 and 20 years.  Supplementary Table 40  reports the association results before and after adjusting for educational attainment. Our analysis shows that the effect sizes shrink after including educational attainment as a covariate, with an average reduction of around 15%. We also estimated the effect of a polygenic risk score of AFB calculated from meta-analysis data excluding the deCODE cohort. The polygenic score remains highly significant. The effect of 1SD of the AFB score decreases from 0.19 years (69 days) without controlling for education to 0.16 years (59 days) when we control for years of education. Second, we estimated how the coefficients change after controlling for Education Attainment (EA) and Age at First Sex using the UKBiobank ( N =50,954). We ran two association models: the first follows the GWAS analysis plan with no additional covariates and the second added years of education and age at first sexual intercourse as covariates. The results are presented in  Supplementary Table 41  and  Supplementary Figure 37 . Our analysis shows that the effect sizes of our top hits are highly concordant (R 2 =0.94). The inclusion of EA and AFS as covariates weakens the effect sizes on average by 40% and increases the  P -value of the estimated coefficients. Overall, we interpret this additional analysis as a robustness test that confirm that the top hits from our meta-analysis are robust to the inclusion of the confounding factors of EA and AFS.\nWe performed a Haploplotter analysis 112  to examine if lead SNPs and/or functional variants identified using RegulomeDB show evidence of positive selection. Three variants showed standardized integrated haplotype scores <-2 or >2, indicating that these variants represent the top 5% of signals in the population. These SNPs are: 1) rs7628058 on chromosome 3 for AFB, an eQTLs for  RBM6  in monocytes; 2) rs2247510 on chromosome 3 for AFB, an eQTL for  RBM6  and  HYAL3  in monocytes and binding site for a range of transcription factors; 3) rs2415984, the lead SNP in the chromosome 14 locus for NEB. Results are presented in  Supplementary Table 42 .","source_license":"Public-Domain","license_restricted":false}