To conduct rare variant burden analyses, we obtained whole exome sequencing (WES) data for 454 787 individuals provided by the UK Biobank (UKBB) study [35]. Individuals were excluded from subsequent analyses if they exhibited excessive heterozygosity, had autosomal variant missingness on genotyping arrays of ≥5%, or were not encompassed within the subset of phased samples as defined by Bycroft et al [36]. Additionally, we excluded individuals who did not have a broadly European genetic ancestry, resulting in a final WES dataset of 420 162 individuals. Additional exclusions occurred during the analysis of specific phenotypes due to incomplete data.

We used the UKBB Research Analysis Platform (RAP; https://ukbiobank.dnanexus.com/) for variant quality control and annotation. The RAP is a cloud-based computing environment and central data repository for UKBB WES and phenotypic data. Customized applets designed for the RAP facilitated additional quality control on exome sequencing data [37], expanding on the procedures outlined in Backman et al [35]. Using provided population-level Variant Call Format (VCF) files, we first split and left-corrected multi-allelic variants into separate alleles using “bcftools norm.” Next, we performed genotype-level filtering using “bcftools filter” separately for single nucleotide variants (SNVs) and insertions/deletions (InDels) using a missingness-based approach. With this approach, SNV genotypes with depth < 7 and genotype quality < 20 or InDel genotypes with a depth < 10 and genotype quality < 20 were set to missing (ie, ./.). We further tested for an expected alternate allele contribution of 50% for heterozygous SNVs using a binomial test; SNV genotypes with a binomial test P value ≤ 1 × 10⁻³ (indicative of a genotype calling error or somatic mutation) were set to missing. Following genotype-level filtering we recalculated the proportion of individuals with a missing genotype for each variant and filtered out all variants with a missingness value > 50%.

Following that, we annotated variants by predicting their functional consequences using ENSEMBL Variant Effect Predictor (VEP) v104 [38]. We applied the “–everything” flag and activating plugins for REVEL [39], CADD [40], and LOFTEE [41]. For each variant, we assigned priority to a specific ENSEMBL transcript based on whether the annotated transcript was protein-coding, MANE select v1.0 [42], or the VEP Canonical transcript, respectively. The consequence for each individual variant was determined by prioritizing severity according to VEP's definitions. After the annotation process, we consolidated stop gained, frameshift, splice acceptor, and splice donor variants into a unified category known as protein truncating variants (PTVs). The consequences for missense and synonymous variants were determined following the definitions provided by VEP. Subsequently, only autosomal or chrX variants within ENSEMBL protein-coding transcripts and within transcripts covered by the UKBB WES assay were retained for subsequent burden testing.

Here, we used the previously reported impact of rare predicted deleterious variants in the NR5A1/SF-1 gene on age at natural menopause (ANM) in this cohort [43], to optimize our definition of deleterious variant class for association testing. We then explored the effect of high-confidence PTVs and missense variants with CADD score ≥25 on other potentially relevant traits in UKBB: age at menarche in women (AAM), age at voice breaking in men (AVB), birth weight, blood pressure (diastolic and systolic), body fat percentage, body mass index (BMI), childlessness, circulating IGF-1 (insulin-like growth factor 1) levels, glucose levels, glycated hemoglobin (HbA1c), high-density lipoprotein (HDL) cholesterol, height, low-density lipoprotein (LDL) cholesterol, number of live births, obesity (BMI >30 kg/m²), red blood cell count, relative adiposity in childhood, SHBG (sex hormone binding globulin) levels, testosterone levels, triglycerides, type 2 diabetes, and waist hip ratio (Supplementary Table S1) [44].

To conduct rare variant burden tests, we used a custom implementation of BOLT-LMM v2.3.6 [45] for the RAP. BOLT-LMM requires 2 main inputs: (i) genotypes at variants with a minor allele count > 100 obtained from genotyping arrays to construct a null linear mixed effects model; and (ii) a larger set of variants collapsed on ENSEMBL transcript to perform association tests. For the former, we extracted genotyping data available on the RAP, limiting it to the same set of individuals used for rare variant association tests. For the latter, since BOLT-LMM expects imputed genotyping data rather than per-gene carrier status, we created dummy genotype files where each variant represented one gene, and individuals with a qualifying variant within that gene were coded as heterozygous, irrespective of the number of variants they carried in that gene.

We generated dummy genotype files for 3 overlapping variant classes: (i) high-confidence PTVs defined by LOFTEE; (ii) missense variants with CADD ≥ 25; and (iii) damaging variants, which included both high-confidence PTVs and missense variants with CADD ≥ 25. For each phenotype tested, BOLT-LMM was executed with default parameters, with the inclusion of the “lmmInfOnly” flag. To obtain association statistics for individual markers, we provided all 26 657 229 markers, regardless of their filtering status, as input to BOLT-LMM. All tested phenotypes were treated as continuous traits or binary traits, adjusted for age, age squared (age²), sex, the first 10 genetic principal components as calculated in Bycroft et al [36], and WES batch as a categorical covariate (50k, 200k, or 450k).

For the phenotype derivation and processing refer to Supplementary Table S1 [44].

To explore how rare deleterious missense variants disrupt SF-1 function, we next categorized variants by protein domain. Briefly, we included predicted deleterious missense variants with CADD score ≥ 25 with MAF < 0.1% which map to either of the 2 main functional domains of SF-1 protein (DNA binding domain [DBD] and ligand binding domain [LBD]), as defined by UniProt. We implemented a generalized linear model (GLM) using the Python package “statsmodels” [46]. For the domain-level burden test, the number of variant alleles across the 2 domains was summed up into a single score under a simple additive model. This score was used as a predictor of the traits in a three-step regression. All traits assessed with BOLT-LMM were included in the GLM burden test, and in addition, we expanded the analysis to encompass additional traits, including spleen volume, platelet count, and reticulocyte count (Supplementary Table S1) [44]. For all association tests we corrected for age, age², the first 10 genetic principal components provided by Bycroft et al [36], and study participant WES batch as a categorical covariate.

All analyses and data manipulations were conducted in R (v4.1.2) and in the UKB Research Analysis Platform (RAP; https://ukbiobank.dnanexus.com/).

Using WES data from 227 858 women of European genetic ancestry, we confirmed that individuals carrying rare deleterious missense variants (with CADD score ≥ 25) in the NR5A1/SF-1 gene had earlier ANM (beta = −1.76 years, [−2.50, −1.01], P = 3.80 × 10⁻⁶, N = 131) (Fig. 1, Supplementary Table S2) [44]. Notably, we did not identify any carrier of rare PTVs alleles in NR5A1/SF-1.

Using the same variant mask of rare deleterious missense variants with CADD score ≥ 25 (Supplementary Table S3) [44] to explore associations with additional traits (Supplementary Table S1) [44], we observed nominal associations of the NR5A1/SF-1 gene with obesity (odds ratio [OR] = 1.041, [1, 1.078], P = .046, N = 472) (Fig. 2), HbA1c levels (beta = .54, [0.11, 0.97], P = .013, N = 451), HDL (beta = .04, [0.012, 0.073], P = 6.5 × 10⁻³, N = 395), circulating IGF-1 levels (beta = −.48, [−0.95, 0], P = .05, N = 415) and systolic blood pressure (beta = .46, [−1.41, 2.32], P = .03, N = 442) (Supplementary Table S2) [44]. Surprisingly, we found no association between this variant mask with other expected reproductive traits, eg, puberty timing in either sex (Supplementary Table S2) [44].

Next, we conducted individual rare deleterious NR5A1/SF-1 variant-level association tests with ANM. We found 7 variants that showed nominally significant individual associations (P ≤ .05) with earlier ANM. Notably, one of these variants (c.43G>A; p.Val15Met, Allele count (AC) = 1) has previously been reported in rare cases of 46,XY sex reversal [47], and premature ovarian insufficiency [48] (Supplementary Table S4) [44]. A second variant (c.721C>T; Arg241Trp, AC = 2) has been also reported in a rare case of 46,XY sex reversal [49], (Supplementary Table S4) [44]. Lastly, a third variant (c.1063G>A; p.Val355Met, AC = 135) has been reported in a boy with bilateral anorchia, micropenis, and progressive testicular degeneration. His mother, who also carried this variant, had 2 spontaneous miscarriages, and had left ovariectomy and removal of the corresponding fallopian tube due to ovarian cysts at the age of 22 (Supplementary Table S4) [44, 50].

To explore how rare deleterious missense NR5A1/SF-1 variants associated with ANM disrupt its protein function, we categorized them by the SF-1 protein domains. Only 3 variants (n = 24 carriers) were found in the hinge region of the SF-1 protein, and none showed an association with ANM (all P > .05). Thus, we focused only on variants in the DBD and LBD.

Using GLM domain-level burden testing, we found that carriers of rare deleterious missense NR5A1/SF-1 variants in each domain had earlier ANM (DBD: beta = −4.03 years, [−5.71, −2.34], P = 2.80 × 10⁻⁶, N = 27; LBD: beta = −1.80, [−2.78, −0.82], P = 3.18 × 10⁻⁴, N = 80). We also tested if considering both DBD and LBD simultaneously increased our strength of association with ANM and found that the association between NR5A1/SF-1 and ANM was strengthened (beta = −2.36, [−3.20, −1.51], P = 4.6 × 10⁻⁸, N = 107). This combined analysis demonstrates a more robust and significant relationship between NR5A1/SF-1 variants and ANM compared to the domain-agnostic test, underscoring the importance of domain-specific effects in influencing NR5A1 function and menopausal timing.

We used this “domains” mask of grouped variants in DBD and LBD to search for NR5A1/SF-1 gene-level associations with other traits beyond ANM, using GLM burden testing. Associations with domain-unrestricted variants are shown for comparison in Supplementary Table S5 [44]. Our analysis revealed that rare deleterious variants in the DBD and LBD of the NR5A1/SF-1 gene were associated with a higher risk of obesity (sex combined: OR = 1.061, [1.01, 1.106], P = .015, N = 344), particularly in women (OR = 1.057, [1.034, 1.163], P = 3.87 × 10⁻³, N = 176), but not in men (OR = 1.019, [0.955, 1.088], P = .57, N = 168) (Supplementary Table S5) [44]. Notably, 32.3% of women who carried rare deleterious missense NR5A1/SF-1 variants located in the DBD and LBD had obesity, compared to 23% of noncarriers. Additionally, we found a positive association of these variants with higher platelet count (beta = 8.55 × 10⁹ cells/L, [2.36 × 10⁹, 14.73 × 10⁹], P = 6.8 × 10⁻³, N = 334) in both sexes (Supplementary Table S5) [44].