Genome-wide analysis identifies 12 loci influencing human reproductive behavior.

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 Until 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.

Genome-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. Cohort 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 . The 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 . We 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 . As 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 ) We 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 ). In 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 . We 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 ). We 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). Following 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 -value0.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 ). We 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). For each of the 12 SNPs identified in the GWAS, local ( cis , 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). We 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 ). Potentially 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 ). DEPICT 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 Previous 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. To 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. We 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, 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 .

We 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 . Associations 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. We 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 ). As 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. To 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 ). To 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 ). We 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). Information 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 . Identifying 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 . The 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. We 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 ). An 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. Five 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. Functional 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. To 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). Several 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. We 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.

This 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. Two 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. The 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 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 Functional 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 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 AFB 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 Our 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.