Understanding the genetic complexity of puberty timing across the allele frequency spectrum

Abstract

Pubertal timing varies considerably and is associated with later health outcomes. We performed multi-ancestry genetic analyses on ~800,000 women, identifying 1,080 signals for age at menarche. Collectively, these explained 11% of trait variance in an independent sample. Women at the top and bottom 1% of polygenic risk exhibited ~11 and ~14-fold higher risks of delayed and precocious puberty, respectively. We identified several genes harboring rare loss-of-function variants in ~200,000 women, including variants in ZNF483, which abolished the impact of polygenic risk. Variant-to-gene mapping approaches and mouse gonadotropin-releasing hormone neuron RNA sequencing implicated 665 genes, including an uncharacterized G-protein-coupled receptor, GPR83, which amplified the signaling of MC3R, a key nutritional sensor. Shared signals with menopause timing at genes involved in DNA damage response suggest that the ovarian reserve might signal centrally to trigger puberty. We also highlight body size-dependent and independent mechanisms that potentially link reproductive timing to later life disease.

Fig. 1. AAM GWAS and gene prioritization.

a, Miami plot showing signals from the European meta-analysis for AAM (top), and G2G scores with names of the top 50 genes annotated (bottom). Top, y axis is capped at −log₁₀(1 × 10⁻¹⁵⁰) for visibility. b, The 50 top-scoring genes implicated by G2G, annotated by their sources of evidence. Relevant data are included in Supplementary Tables 2 and 13–16.

Fig. 2. Exome-wide rare variant burden associations with AAM.

a, Manhattan plot showing gene burden test results for AAM. Six genes passing exome-wide significance (P < 1.54 × 10⁻⁶) are highlighted. Point shapes indicate variant-predicted functional class (DMG; damaging). b, Quantile–quantile (QQ) plot for gene burden tests. c, Comparison of gene burden associations (effect sizes with 95% CIs) for AAM (female participants, in years, n = 222,283) and AVB (male participants, in three categories, n = 178,625). Relevant data are included in Supplementary Table 5.

Fig. 3. Epistatic interactions between rare coding variants and common genetic susceptibility on AAM in the UK Biobank.

a, Interaction effects (mean and 95% CIs) on AAM between a GWAS PGS and carriage of qualifying rare variants in the six exome-highlighted genes (n = 222,283). b,c, Predicted mean (with 95% CI) AAM in (b) noncarriers (black) and carriers (light blue) of rare variants in six genes without significant (P > 0.05) interaction effects and (c) in noncarriers (left) and carriers (right) of rare variants in ZNF483, which shows significant interaction (P = 0.025). In c, points show individual AAM values, with the color gradient indicating increasing age. d, Plot of individual rare damaging (DMG) variant associations with AAM by ZNF483 functional domains. The coding part of ZNF483 is depicted by the horizontal black line. Included DMG variants had an MAF < 0.1% and were annotated to be either HC PTVs or missense variants with CADD score ≥25. Each variant association is represented by a circle and vertical line—the line length indicates the −log₁₀(P), in the direction of its effect on menarche in carriers of the rare allele, and the circle size indicates the number of carriers of each variant (that is, allele count). Relevant data are included in Supplementary Table 12.

Fig. 4. Stratification of AAM signals and biological pathway enrichment by their influence on early childhood weight.

a, Proportion of GWAS signals for AAM by early childhood weight trajectory. b, Biological pathways enriched for high-confidence AAM genes (left), plus enrichment within early childhood weight trajectories (right). Pathway enrichment was calculated using g:Profiler, and the displayed P values are Bonferroni-corrected. Row names describe pathway clusters. Strength of associations with individual pathways is indicated by circles. Circle size reflects the proportion of pathway genes that are high-confidence AAM genes (in percentage). Right, determining whether each pathway cluster remains enriched for AAM genes when stratified by early childhood weight trajectory. Supporting data are included in Supplementary Tables 21 and 23–26.

Fig. 5. Enrichment of gene drivers of GnRH migration and maturation in the AAM GWAS.

a, Schematic representation of the stages of GnRH neuron migration during embryonic development. Using RNA-seq data, Zouaghi et al. grouped differentially expressed genes into 23 expressional trajectories based on their comparative level of expression during the early (yellow), intermediate (amber) and late (red) stages of GnRH migration. b, Genome-wide MAGMA enrichment for AAM gene associations within each expression trajectory. Colored symbols accompanying the points indicate the expression level of genes within each trajectory at the three aforementioned stages of migration; that is, the first yellow symbols indicate whether gene expression remains the same (=), is upregulated (^) or downregulated (v) in the early stage of migration. Orange and red symbols denote the same for the intermediate and late stages of migration, accordingly. c, Trajectories significantly (P_adj < 0.05) enriched at the genome-wide level in b are highlighted in red and also show significant (P < 0.05) overlap with the 665 high-confidence AAM genes. Supporting data are included in Supplementary Tables 18 and 27.

Fig. 6. Interactions between GPCRs on AAM.

a, A total of 24 brain-expressed GPCRs were implicated in AAM by the G2G analysis of the white European GWAS data. Point colors indicate the number of concordant G2G predictors implicating each GPCR. b, Time-resolved NDP-αMSH-stimulated cAMP production in HEK293 cells expressing MC3R-alone or with both MC3R and GPR83. Data are presented as the mean (±s.e.) percentage of the maximal MC3R-alone response (from six independent experiments). c, AAM according to dosage of MC3R function-increasing C alleles at rs3746619 (x axis in each panel) and GPR83 expression-increasing T alleles at rs592068. Predicted means are represented by the lines, and the accompanying 95% CIs are denoted by the shaded areas. β_interaction = −0.034 ± 0.015 years, P_interaction = 0.02. Supporting data are included in Supplementary Tables 28 and 30.