Identification of novel variants underlying non-syndromic primary ovarian insufficiency using a targeted NGS gene panel.

Primary ovarian insufficiency (POI) is a complex disease affecting approximately 1-4% of women worldwide ( 1 ). POI is characterized by amenorrhea in women under of the age of 40, hypoestrogenism, and elevated levels of follicle-stimulating hormone (FSH), leading to progressive loss of ovarian function and infertility ( 2 – 4 ) The long-term consequences of POI include an increased risk of osteoporosis, cardiovascular disease, and mental health disorders ( 5 , 6 ). As such, POI has a substantial impact on women’s reproductive and overall health. The etiology of POI is heterogeneous and may result from various causes, including iatrogenic and environmental exposures ( 7 , 8 ) autoimmune and metabolic diseases ( 9 , 10 ), viral infection ( 11 ), X-linked abnormalities ( 12 , 13 ) and autosomal gene mutations ( 14 , 15 ). To date, the genetic mechanism underlying POI remain poorly understood, no standardized diagnostic methods currently exist to detect its genetic basis. POI may present as part of a pleiotropic genetic syndrome or as an isolated, non-syndromic form. The advent of the “omics” era has led to the discovery of several genes implicated in both syndromic and non-syndromic POI due to their key roles in reproductive function ( 16 – 28 ). Given the high genetic heterogeneity and relatively low prevalence of POI, traditional approaches such as candidate gene studies and genome-wide association studies (GWAS) are limited in their ability to elucidate the underlying genetic causes. While some pathogenic variants have been identified in specific families or subsets of patients, their contribution to the POI phenotype often remains uncertain. This knowledge gap delays the development of accurate diagnostic tools and personalized treatment strategies. To address these limitations, next-generation sequencing (NGS) has emerged as a powerful technology for the comprehensive analysis of genes involved in ovarian function ( 21 ). Recent studies ( 29 – 32 ) have demonstrated the utility of targeted NGS panels for identifying potentially pathogenic variants, allowing the simultaneous analysis of multiple genes and improving diagnostic efficiency in terms of cost and turnaround time. Targeted NGS represents a promising strategy to investigate selected genes associated with ovarian biology and enhance clinical management through more precise etiologic diagnoses ( 33 ). In this study, we analyzed the genetic profile of women with non-syndromic POI using a targeted Ampliseq panel targeting the exons of 72 genes, selected based on a previous systematic review conducted by our group ( 34 ). The aim was to investigate the genetic basis of POI in a cohort of Italian patients selected according to OvAge © algorithm ( 35 ). The OvAge algorithm integrates clinical, hormonal, and ultrasound parameters into a mathematical model that yields a single value, expressed as “ovarian age.” Its main advantage is the combination of complementary markers, which reduces the variability and limitations of individual tests and provides an intuitive output that can be directly compared with chronological age. Our ultimate goal was to develop a diagnostic test that integrates genetic information for clinical use, supporting the early identification of women at risk of POI and enabling timely fertility preservation strategies.

The aim of this study was to characterize the mutational profile of individuals affected by non-syndromic primary ovarian insufficiency (hereafter referred as nsPOI). To this end, we sequenced all coding exons of 72 genes using the Ampliseq technology on the Ion Torrent platform ( Supplementary File S1 ). The experimental workflow is outlined in Figure 1 . Workflow of the NGS analysis. The flowchart shows the workflow to identify the potential causative genetic variants in our nsPOI patients. Legend: AF, allele frequency; P, pathogenic; LP, likely pathogenic; VUS, variant of uncertain significance; LB, likely benign. Genomic DNA from peripheral blood samples of 100 Italian women with nsPOI and 200 healthy controls was sequenced (see Materials & Methods for inclusion and exclusion criteria). Variants were filtered to exclude those with a frequency greater than 1% or those not present in the non-Finnish European population of the gnomAD and 1000 Genomes databases ( Figure 1 ). By comparing the nsPOI cohort with healthy controls ( Supplementary File S4 ), we identified 78 unique variants ( Supplementary File S5 ) that were exclusively present in 60 of the 100 nsPOI patients. The remaining 40 patients did not carry any mutation in the analyzed genes ( Figure 1 ). These variants included 70 missense (89,7%), 1 nonsense (1,3%), 1 frameshift (1,3%), 5 in-frame deletions (6,4%) and 1 in-frame insertion (1,3%) variants ( Figure 2A ). These 78 variants were distributed across in 41 genes ( Supplementary File S5 ) encoding 14 transcription factors, 11 proteins involved in signal transduction and the cell cycle, 10 meiotic factors and 6 enzymes ( Figure 2B ). Genes involved in meiosis showed the highest mutation frequency (31%), followed by those involved in signal transduction and cell cycle (26%), transcriptional regulation (23%) and enzymatic activity (20%). The Gene Ontology analysis of the mutated genes is presented in Figure 2C and Supplementary File S6 . The most enriched functional categories (those including >6 genes) were related to the development of the female reproductive system and gamete generation, although many genes were also involved in transcriptional regulation and maintenance of the stem cell niche. Genomic characteristics of the nsPOI cohort. (A) Pie chart showing the frequency of each mutation type identified in nsPOI patients. (B) Histogram displaying the prevalence of variants per gene. The 78 filtered variants were found in 41 genes, encoding (from left to right): 10 meiotic factors, 11 proteins related to signal transduction/cell cycle, 14 transcription factors, and 6 enzymes. The y-axis indicates the frequency of variants in each gene (blue) and the percentage of patients carrying variants in that gene (pink). Percentages of each gene class correspond to the sum of blue bars. (C) Tree plot showing hierarchical clustering of positively enriched GO terms among the 41 mutated genes. Circle size reflects the number of genes per GO term; circle color indicates adjusted p-value significance. Clusters are annotated with representative keywords. (D) Pie chart presenting the proportion of patients with LP/P or VUS variants versus those with no detectable mutations (negative). To classify the identified variants, we followed the ACMG guidelines ( 41 ), using five predictive algorithms (SIFT, PolyPhen-2, MutationTaster, CADD, and Grantham) for optimal interpretive accuracy (see Materials and Methods and Supplementary File S3 for details). Variants were classified as likely benign/benign (LB/B), likely pathogenic/pathogenic (LP/P) and of uncertain significance (VUS). In total, we identified 2 LB variants, 24 LP variants, 2 P variants, 50 VUS variants. Overall, 40% of patients were mutation-negative, while 60% carried at least one variant. Among these, 23% had LP/P variants and 37% had only VUS ( Figure 2D ). Representative Sanger sequencing validations of selected mutations are shown in Supplementary Figure S1 . The protein Dynein axonemal heavy chain 5 encoded by the DNAH5 gene was mutated in twelve patients (12%). Three patients carried the S3774P variant, two had I3568T, two patients carried the variant T806I. Seven additional patients each presented a distinct mutation: A769V, S2605L, N1420D, L1339R, N934S, Q2949E, and Y4308C. Notably, S3774P and I3568T co-occurred in two patients. Except for T806I, N934S, L1339R, and I3568T, which were located in linker domains, the remaining mutations resided in conserved functional domains ( Figure 3A ). Mutations A769V and N1420D were located in the N-terminal region 1 (DHC_N1) and region 2 (DHC_N1) of the dynein heavy chain, respectively; mutations S2605L and Q2949E were located in the P-loop containing the dynein motor region D3 (AAA_7) and D4 (AAA_7), respectively; mutation S3774P was located in the ATP-binding dynein motor region D5 (AAA_9); and Y4308C was located in the dynein heavy chain region D6 of the dynein motor ( Figure 3A ). Mutations S2605L and N934S were classified as likely pathogenic, while the remaining mutations were categorized as variants of uncertain significance (VUS) ( Table 1 ). Mutation mapper plots of the most frequently mutated genes in nsPOI patients. The plots show the distribution of variants in DNAH5 (A) , LAMC1 (B) , HK3 (C) , HSD17B4 (D) , and PCDH11X (E) . Amino acid changes are indicated above each lollipop, and line length represents the number of patients carrying the corresponding variant. Plots were generated using the Mutation Mapper tool on the cBioPortal platform ( https://www.cbioportal.org/mutation_mapper ). Variants identified by targeted NGS in our POI cohort. The table shows the most frequent variants identified in 60/100 patients in our cohort. For each variant the allele frequency in population databases (GnomAD, 1000 genomes_EUR), the pathogenicity prediction scores (Sift, Polyphen, Mutation Taster, CADD, Grantham) and the final classification according to ACMG criteria are reported. Abbreviations: GnomA NEF, Genome Aggregation Database Non-Finnish European; Eur AF, European population Allele Frequency; B, benign; D, deleterious; NA, not available; P, pathogenic; LP, likely pathogenic; VUS, variant of uncertain significance; LB, likely benign. Six patients had 4 missense mutations (R1011H, Y1035S, A1239V, A1335S) in the protein encoded by the laminin subunit gamma 1 ( LAMC1 ) gene. The R1011H mutation, which is localized in the laminin EGF domain ( Figure 3B ), was present in 3 patients. All mutations had a very low frequency or were not present in the public control databases, and their prediction scores supported potential pathogenicity at different levels ( Table 1 ). The gene encoding hexokinase 3 ( HK3 ) was mutated in five patients carrying 4 missense mutations of the gene (C237R, I347T, Q600H, P676S). The mutations C237R and I347T are located in the hexokinase 2 domain, while the mutation Q600H is located in the hexokinase 1 domain and the mutation P676S is located in the junction domain between the hexokinase 1 and hexokinase 2 domains at the C-terminus ( Figure 3C ). In particular, C237R and I347T are likely pathogenic variants ( Table 1 ). Four patients had four different missense mutations of the HSD17B4 gene encoding 17-beta-hydroxysteroid dehydrogenase 4, which is involved in the peroxisomal beta-oxidation pathway for fatty acids. The I53M and C214S mutations were located in non-functional domains, while the R658H and A741S missense mutations were located in the sterol carrier protein 2 (SCP2) domain, which is responsible for the binding and transfer of sterols and phospholipids ( Figure 3D ). The mutations I53M and A741S were classified as likely pathogenic variants ( Table 1 ). In our cohort, four patients had missense mutations in PCDH11X. The mutation T790S, located in the protocadherin domain, was found in three patients and R1010I, at the C-terminus of the protein, was present in one patient and was predicted likely pathogenic ( Figure 3E , Table 1 ). We also identified mutations in known POI-related genes involved in folliculogenesis and ovarian function, many of which have been validated in mouse models ( 49 ). Several X-linked genes are known contributors to POI pathogenesis ( 50 , 51 ). In addition to PCDH11X, we identified Androgen receptor (AR) mutations in three patients. Two distinct deletions (Q66–80del and Q74–80del) within the poly-Q repeat at the N-terminus were classified as pathogenic ( Figure 4A , Supplementary Figure S1 ). Mutation mapper plots of genes with LP/P variants in nsPOI patients. The plots show the distribution of mutations in AR (A) , ADAMTS19 (B) , FANCM (C) , FANCC (D) , GDF9 (E) , SPIDR (F) , and SF1 (G) . Amino acid changes are indicated above each lollipop, and line length reflects the number of patients carrying the corresponding variant. Plots were generated using the Mutation Mapper tool on the cBioPortal platform ( https://www.cbioportal.org/mutation_mapper ). Mutations were also found in ADAMTS1 and ADAMTS19, which encode metalloproteinases involved in ovarian extracellular matrix remodeling ( 52 ). Three patients carried R64C, L117V, and G202S missense variants in ADAMTS19 ( Figure 4B ). The R64C mutation, located in the pro-domain crucial for enzyme folding and activity, was deemed likely pathogenic ( Table 1 ). Additionally, six patients carried VUS in ADAMTS1, further discussed below. In recent years, the Fanconi anemia complementation group ( FANC ) genes have also been implicated in POI ( 22 , 23 ). Mutations were identified in FANCM (3%), FANCC (2%), and FANCG (1%) in our cohort. In FANCM, three missense mutations (A48D, L57F, P1255L) were found ( Figure 4C , Table 1 ). Although not located in known domains, A48D and L57F were predicted likely pathogenic. FANCC mutations included a potentially pathogenic E273Q substitution and a frameshift mutation (I121Tfs*7) causing a premature stop codon ( Figure 4D , Table 1 ). One patient carried the W122C mutation in the FANCG protein ( Table 2 ). This substitution occurred at an interspecies conserved residue, which given the physicochemical distance between the W and C residues (Grantham score=215) should result in a severe disruption of the protein structure, although it was classified as VUS. Mutations in genes carrying only VUS. FOR each variant the allele frequency in population databases (GnomAD, 1000 genomes_EUR), the pathogenicity prediction scores (Sift, Polyphen, Mutation Taster, CADD, Grantham) and the final classification according to ACMG criteria are reported. Abbreviations: GnomA NEF, Genome Aggregation Database non-Finnish European; Eur AF, European population Allele Frequency; B, benign; D, deleterious; NA, not available; VUS, variant of uncertain significance. The gene encoding growth differentiation factor 9 ( GDF9 ), crucial for folliculogenesis and oocyte development, was mutated (R454C) in two patients. The R454C mutation occurs at a conserved amino acid position and is likely to have functionally detrimental effects ( Figure 4E , Table 1 ). Two patients carried mutations in the Scaffolding Protein Involved in DNA Repair ( SPIDR ) gene, which encodes a scaffold protein involved in HR repair. Both a nonsense mutation (S64*) and a missense mutation (R294K) occurred in the DUF4502 domain, the role of which is not yet fully understood ( Figure 4F , Table 1 ). The mutation S64* introduces a premature stop codon near the N-terminus and was classified as highly pathogenic. Two mutations were found in the splicing factor 1 (SF1) gene ( Figure 4G , Table 1 ). A deletion G39_P47del, in the helix-hairpin domain at the N-terminus of SF1 ( 53 ) was classified as LP. This helix-hairpin domain is required for cooperative recognition of 3’ splice sites by stabilizing a unique quaternary arrangement of the SF1-U2AF65 RNA complex during assembly of the spliceosome ( 54 ). In contrast, G517A is categorized as VUS in the low complexity region. Finally, we identified likely pathogenic variants of SOHLH1, AMHR2, TP63, NANOS3, MRPS22, FOXL2 and DMC1 genes in one patient each. These mutations were absent or novel in public databases and had consistently high pathogenicity scores across prediction tools. Of particular note, we identified a novel Gly323 deletion in SOHLH1, within a low-complexity region involved in folliculogenesis, where other deleterious missense variants have been described ( 55 ). In addition, one patient carries a mutation of the TP63 gene, which codes for a transcription factor of the p53 family. TP63 is expressed in primordial and primary follicles and preserve the germ line integrity ( 56 ). The R487C mutation is located in the linker domain but predicted as potentially deleterious ( 57 , 58 ). Our analysis also identified additional 20 genes that carried only VUS. As previously described, six patients harbored three missense mutations (T514A, T732I, A806V) in another member of the ADAMTS family, ADAMTS1 ( Table 2 ). Two patients carried the T514A mutation, located in a linker domain ( Figure 5A ). The T732I variants, present in 3 patients, and the A806V were both located within the ADAM spacer domain ( Figure 5A ). Notably, T732I had the highest Grantham score among the three, indicating that the substitution of threonine—a small, polar, hydrophilic amino acid—with isoleucine—a large, non-polar, hydrophobic residue—could significantly alter protein structure. Mutation mapper plots of genes with VUS variants in nsPOI patients. The plots show the distribution of mutations in ADAMTS1 (A) , BNC1 (B) , CDKN1B (C) , WT1 (D) , and WDR62 (E) . Amino acid changes are indicated above each lollipop, and line length corresponds to the number of patients carrying the respective variant. Plots were generated using the Mutation Mapper tool on the cBioPortal platform ( https://www.cbioportal.org/mutation_mapper ). CDKN1B and WDR62 plots were modified with Biorender . Three patients carried the same variant (R315C) in the POF1B gene, which has been previously associated with POI ( 59 – 61 ). Although classified as VUS, this variant was predicted to be deleterious by four out of five in silico algorithms ( Table 2 ). The BNC1 gene, encoding the zinc finger protein Basonuclin 1, is highly expressed in ovarian germ cells and plays a role in transcriptional regulation ( 62 , 63 ). We identified two missense mutations (E209K and G661R) in the linker domain of the BNC1 protein ( Figure 5B , Table 2 ). The E209K mutation was predicted to be deleterious by the SIFT, PolyPhen, CADD and MutationTaster algorithms, likely due to the change in charge affecting structural stability and DNA binding capacity. In the CDKN1B gene, which encodes the cell cycle inhibitor p27^Kip1, two missense variants were identified: S7C, located in the Neurogenin-2 binding domain at the N-terminus ( 64 ), and P117S which replaces a highly conserved proline within the Jab1/CSN5 binding domain ( Figure 5C ). In the WT1 gene, which encodes the WT1 transcription factor, we found two missense mutations ( Figure 5D ). The G37S mutation was located at the N-terminus, which is characterized by a proline- and glutamine-rich DNA-binding domain, while the S325L mutation was located in one of the four zinc finger motifs at the C-terminus (ZnF_CH2H2 domain). Several variants have been described for this gene, including G37S and S325L, both of which are listed with very low frequency in the gnomAD database. Two patients carried two different mutations, S10L and L35I, in the N-terminus of the protein encoded by the MutS homolog 4 (MSH4) gene, which is involved in the mismatch repair process. We found two mutations in the gene encoding the WD repeat-containing protein 62 ( Figure 5E ). The mutations S275L and Y336C were located in the first WD40 domain of MABP1/WDR62 ( 65 ) of the protein, which is involved in the interaction with Aurora A and binding to microtubules at the spindle pole ( 66 ). In addition, we found VUS mutations in the genes SYCE1, AGTR2, FOXO1, FOXO3, BBS9, DUSP22, LARS2, NUPR1, POLR3H, SALL4 in one patient each. See Table 2 and Supplementary File S5 . Among the 60 patients with genetic variants, we found that 37 had alterations in a single gene (monogenic patients), while 23 carried variants in two or more genes (polygenic patients) ( Figure 6A , Supplementary File S7 ). Of the 37 monogenic patients, 27 (73%) carried VUS, while 10 (27%) had likely pathogenic mutations ( Table 3 ). One patient carried the truncating pathogenic variant, which occurs in the SPIDR gene. Co-occurrence of mutations in the nsPOI cohort. (A) The graph shows the percentage of POI patients in whom one or more genes are simultaneously mutated. (B) Stacked bar plot reports the number of mutated genes (y-axis) per patients (x-axis). Gene names are displayed in the bars. The table lists the 37 patients that could be explained by monogenic variants. P, pathogenic; LP, likely pathogenic; VUS, variant of uncertain significance. Genes exclusively mutated in monogenic patients included HK3, NOTCH2, MSH4, BNC1, POU5F1, DUSP22, FANCG, MRPS22, and SYCE1 ( Table 3 ). Four of the 12 patients with DNAH5 mutations had only a VUS mutation of the gene. Among the polygenic cases, 15 patients had mutations in two genes (15%), 6 patients had mutations in three genes (6%) and only 2 patients had mutations in four genes (2%) ( Figure 6B , Supplementary File S7 ). Of particular interest, two patients shared the same combination of mutations in the ADAMTS1 (T732I) and GDF9 (R454C) proteins, while two other patients had concurrent but different mutations in the ADAMTS1 and DNAH5 proteins: ADAMTS1-A806V and DNAH5-S2605L; ADAMTS1-T514A and DNAH5-T806I. All POI patients were also screened for germline CNVs using amplicon-based NGS data generated with Ion Reporter software. A pooled DNA sample from patients in the control group served as reference baseline. Based on the MAPD values, 99 out of 100 patients were eligible for CNV analysis (see Materials and Methods). CNVs were detected in 11% of patients, involving SYCE1, DUSP22, and INHBB genes ( Table 4 ). Chromosome regions with CNVs identified in POI cohort. The table shows CNVs and co-occurred mutations in genes previously associated to POI. For each CNV chromosome region, type of CNV, genomic coordinates CNV length, gene, CADD SV Phred, StrVCTVRE, X-CNV(MVP), co-occurred mutated genes and their amino acid changes scores are indicated. Five patients had CN amplifications of the SYCE1 gene on chromosome 10q26.3 and all were identified as focal (<3Mb). Four of the five patients with copy number gains also had missense mutations and deletions in other genes ( Table 4 ). Notably, one patient with SYCE1 CN gains also had LP variants in AR (Q66_Q80del) and FOXL2 (A234del) ( Table 4 ), while the other had LP mutations of LAMC1 (A1335S) and AR (Q74_Q80del) and a VUS (A140S) in the FOXO3 gene. Five patients had focal CNVs of the DUSP22 gene on chromosome 6p25. Four patients had CN gains of the gene, while one patient had a copy number loss of the gene. Three of the four patients with copy number gain also had a pathogenic mutation each in DNAH5 (Q2949E), LAMC1 (R1011H) and HK3 (C237R) ( Table 3 ). We also found a copy number loss of the INHBB gene on chromosome 2q14.2 in a patient who had no mutations in other genes ( Table 3 ). For pathogenetic classification of the identified CNV, we used the StrVCTVRE, MVP and CADD_SV_Phred scores (see Material and Methods for more details). According to the scores obtained, the microdeletion found in the INHBB gene was classified as likely pathogenic. Supplementary Figure S2 shows representative Q-PCR analysis of CNVs in SYCE1, DUSP22 and INHBB genes.

In this study, we analyzed a cohort of 100 Italian women affected by non-syndromic primary ovarian insufficiency (nsPOI), selected using the OvAge © algorithm ( 35 ), to investigate the contribution of genetic alterations using a targeted next-generation sequencing (NGS) approach. Preliminary validation has shown that OvAge discriminates pathological conditions such as POI and Polycystic Ovary Syndrome (PCOS), which are often misclassified by individual tests ( 35 ). In longitudinal observations, OvAge was also able to anticipate menopause within 1–2 years in women predicted to have an ovarian age close to 50 years, suggesting potential prognostic value. However, it was derived from a single, albeit large, population and requires advanced ultrasound technology, so external validation, broader clinical application, and integration with genetic information are warranted. In the present study, by integrating sequencing data with a panel of 72 POI-related genes, we identified 78 rare variants distributed across 41 genes, which were absent in 200 matched controls. Variants were classified following ACMG guidelines, resulting in a diagnostic yield of 60%, with 23% of patients carrying pathogenic or likely pathogenic (LP/P) variants and 37% harboring only variants of uncertain significance (VUS). Our findings further support the notion that nsPOI is genetically heterogeneous disorder, with both monogenic and polygenic contributions. Notably, DNAH5 emerged as the most frequently mutated gene (12%), with several variants located in conserved motor domains. DNAH5 encodes a microtubule-associated motor protein involved in ciliary motility. Mutations in DNAH5 have previously been linked to primary ciliary dyskinesia (PCD) ( 67 ), as well as non-syndromic asthenozoospermia and hypospadias ( 68 , 69 ). Aboura and colleagues identified chromosomal amplifications at 5p14.3, the locus of DNAH5 , in a POI cohort. To our knowledge, this is the first report describing DNAH5 mutations in nsPOI patients. We also identified mutations in the gene encoding the laminin subunit gamma-1 ( LAMC1 ). LAMC1 is known to interact with other laminin family proteins in the extracellular matrix to promote ovarian follicle development ( 70 , 71 ). Overexpression of LAMC1 is involved in the progression of gynecologic cancers ( 72 , 73 ) and predicts poor prognosis in gastric and esophageal cancers ( 74 – 77 ). In contrast, only one missense mutation ( 78 ) and one SNP haplotype ( 79 ) have been associated with increased POI risk without functional characterization. In this work, we identified four novel missense mutations (R1011H, Y1035S, A1239V, A1335S) in LAMC1 that are very rare or absent in the general population. In particular, the mutation R1011H, which is localized in the EGF_like domain and occurred in three different patients, was classified as likely pathogenic. Our data also highlight the relevance of the ADAMTS family in ovarian function ( 80 , 81 ). We detected pathogenic or VUS variants in both ADAMTS1 and ADAMTS19, proteases known to regulate folliculogenesis and ovulation via extracellular matrix remodeling ( 82 – 84 ). The role of ADAMTS1 and ADAMTS19 in reproductive function has been demonstrated using mouse models: ADAMTS1 homozygous knock-out mice had fewer ovarian follicles ( 85 ), while the protease ADAMTS19 is overexpressed in the gonads of female mice ( 86 ). Regarding the role of ADAMTS19 in nsPOI, so far only SNPs in intronic regions of ADAMTS19 have been potentially associated with POI in both Caucasian ( 87 ) and Asian populations ( 79 ). In our study, we described three novel mutations in ADAMTS1 (A806V, T732I and T514A) found in six patients. In particular, the mutation T732I, which occurred in three patients, is located in the spacer domain, which is necessary for association with ECM components and regulation of enzyme activity ( 88 ). In addition, we identified three missense mutations (R64C, L117V, G202S) of ADAMTS1 9, of which the R64C mutation was classified as likely pathogenic according to the ACMG criteria. The HSD17B4 gene, also known as D-bifunctional protein (DBP), is a bifunctional enzyme involved in the conversion of androstenedione to testosterone and estrone to estradiol. In this study, we found four different missense mutations (I53M, C214S, R658H and A741S) in HSD17B4 gene in 4 different patients. Studies on Perrault syndrome ( 89 ) provided evidence of the importance of HSD174B for healthy ovarian function. In addition, Puyan et al. found a haplotype and two missense SNPs in the HSD17B4 gene associated with susceptibility to POF in a genetic case-control association study ( 79 ). Remarkably, 31% of mutated genes in our study were involved in DNA repair and meiosis, pathways essential for oocyte integrity. We found mutations in three different genes of the FANC family, FANCM, FANCC and FANCG, whose role in repair during HR in meiosis has been recently investigated ( 90 ). Our results support the observations that defects in the FANC genes can impair normal oogenesis. In our cohort, 3 of the five mutations identified in the FANCM and FANCC genes were classified as LP or P ( Table 2 ). In addition, we found mutations in the gene MSH4, which is required for optimal reciprocal recombination and appropriate segregation of homologous chromosomes during meiosis I ( 91 ). Mutations of this gene have recently been associated with failure of gametogenesis in both sexes ( 92 ). A causal role may also be attributed to mutations in the homologous recombination repair gene SPIDR gene, whose alteration has been associated with gonadal disgenesia ( 93 ) and ovarian failure ( 23 ). We have also identified novel variants in genes previously associated with POI such as PCDH11X, AR, TP63, BNC1, WT1, CDKN1B, SYCE1 as well as previously described mutations such as R454C in GDF9, which is considered one of the causative alterations in POI ( 94 ). It is noteworthy that, to our knowledge, no POI-associated missense mutations for PCDH11X have been reported to date, only CNV ( 52 ). In addition, three patients had two different deletions in the poly-Q region of the AR gene. Previous reports have described ( 95 – 97 ) only point mutations of the gene encoding the AR in POI patients, while short poly-Q polymorphisms have been associated with poorer prognosis only in endometrial cancer ( 98 ). Recent data have identified TP63 mutations in syndromic and non-syndromic POI, associated with oocyte apoptosis and early ovarian reserve depletion ( 17 , 61 , 99 ). The R487C mutation identified in our study has previously been linked to colorectal cancer ( 57 ), low-grade gliomas ( 58 ), and more recently, to pyroptosis-related gene networks ( 58 , 100 , 101 ). This supports a possible role for p63 as a high-risk biomarker in cancer and reproductive disorders. Interestingly, the P117S mutation in the cell cycle inhibitor CDKN1B has been associated with multiple endocrine neoplasia type IV (MEN4) ( 102 , 103 ). This mutation is located in the domain of binding to Jab1/CSN5, which promotes the translocation of p27 from the nucleus to the cytoplasm, thereby favoring cell proliferation ( 104 ). POI-associated genes include the Synaptonemal Complex Central Element 1 (SYCE1) gene which encodes a member of the synaptonemal complex that links homologous chromosomes during prophase I of meiosis. Allelic variants of this gene have been associated with premature ovarian failure ( 105 ) and spermatogenic failure. We found a VUS in this gene (R59K) that was reported as a somatic mutation in a patient with malignant skin cancer ( 106 ). Of particular interest, several variants were shared among multiple patients, and distinct mutations in the same gene were observed across the cohort. In 37 patients, POI could be attributed to monogenic variants, while in 23, mutations in two or more genes suggested a polygenic etiology. We found that two patients had a specific combination of mutations in the ADAMTS1 (T732I) and GDF9 (R454C) genes. Two patients carried both DNAH5 variants, S3774P and I3568T. As parental DNA was unavailable, phasing could not be established, and compound heterozygosity for an autosomal recessive mechanism remains unconfirmed. Segregation studies in relatives and analysis of a larger cohort are planned to assess the contribution of this co-mutation to nsPOI. CNV analysis identified alterations in SYCE1, DUSP22 , and INHBB . Of these, only the INHBB CNV was classified as likely pathogenic by prediction tools. Although SNV mutations in the INHBB gene were reported in a previous study ( 107 ), no evidence of POI-causing CNV in this gene has been described to date. INHBB encodes the βB subunit of inhibin/activin dimers which regulate granulosa cell proliferation and folliculogenesis ( 108 , 109 ). Loss of INHBB function may impair activin-mediated signaling, leading to defective follicular maturation and contributing to premature ovarian failure. We speculate that CNVs involving DUSP22 may dysregulate JNK signaling ( 110 ), thereby promoting apoptosis during follicle maturation. Collectively, these findings suggest that CNVs may play a significant role in the pathogenesis of nsPOI. From a clinical perspective, identification of pathogenic variants—including both SNVs and CNVs—offers opportunities for translation into patient care. Preimplantation genetic testing (PGT) could help prevent transmission of deleterious variants in families with a history of POI. Moreover, women at increased risk could be counseled regarding early fertility preservation strategies, such as oocyte cryopreservation. As the disrupted signaling pathways become more clearly defined, these insights may also pave the way for targeted therapeutic approaches, including pharmacological modulation and, in the future, gene-based interventions. In summary, we employed a customized targeted NGS panel to analyze DNA from 100 POI patients and 200 controls. The approach yielded a diagnostic rate of 60%, implicating 42 genes. These findings provide new insights into the complex genetic architecture of POI. We acknowledge that the pathogenic potential of many variants identified requires confirmation through functional validation in cell and animal models. Several VUS in our cohort appear to be strong candidates for reclassification as “likely pathogenic.” Consistent with recent reports, a substantial proportion of VUS are ultimately reclassified, although the average time to reclassification is approximately 2–3 years, underscoring the need for periodic reanalysis ( 111 – 113 ). Notably, variants such as BNC1 E209K, CDKN1B P117S, WT1 S325L, WDR62 S275L, FANCG W122C, SOHLH2 S147L, and SYCE R59K showed high deleteriousness scores, were located in conserved domains, and were absent or extremely rare in population databases. However, additional functional or segregation evidence will be required to support their reinterpretation. Expanding the cohort size may further improve sensitivity and strengthen genetic associations. We are acknowledge the absence of segregation testing in family members, which could have clarified the significance of several VUS. This was primarily due to budgetary constraints, as the project was supported by a ministerial grant that covered only affected patients and did not extend to relatives. Despite these limitations, the integration of clinical, hormonal, and ultrasonographic data (via the OvAge algorithm) with targeted NGS represents a promising strategy for early POI diagnosis and personalized management. While available treatments can assist with fertility, the irreversible depletion of ovarian reserve in POI remains incurable. As affected women are also at increased risk for comorbidities that reduce life expectancy, early identification of genetic risk could improve long-term outcomes and quality of life.

A cohort of 100 women with non-syndromic POI was recruited from the OvAge database, maintained at the Department of Experimental and Clinical Medicine of the Magna Graecia University of Catanzaro. This database contains data on more than 1,000 women, including serum levels of anti-Müllerian hormone (AMH), follicle-stimulating hormone (FSH), estradiol (E2), three-dimensional antral follicle count (AFC), vascular index (VI), flow index (FI), and vascular flow index (VFI). For all patients, ovarian reserve had previously been estimated using our proprietary OvAge algorithm ( 35 ). This algorithm is a mathematical formula that combines biochemical and ultrasonographic parameters to produce a single interpretable value, referred to as “OvAge”—an estimate of ovarian age—calculated using the following linear equation: OvAge = 48.05 – 3.14 × AMH + 0.07 × FSH – 0.77 × AFC – 0.11 × FI + 0.25 × VI + 0.1 × AMH × AFC + 0.02 × FSH × AFC. The formula was generated through the use of a Generalized Linear Model (GzLM). For this study, patients were selected from the database based on the following criteria: onset of oligomenorrhea or menopause before the age of 40, elevated FSH level (> 40 mU/L), estradiol (E2) < 20 pg/ml, anti-mullerian hormone (AMH) < 1 ng/ml. Additionaly, patients were included if the difference between their OvAge and chronological age exceeded 10 years. A control group of 200 healthy women with an OvAge–chronological age difference of ±2 years was also included. The following exclusion criteria were applied to both POI patients and controls: use of estrogens or progestins or breastfeeding within two months prior to enrollment, ongoing pregnancy, history of endometriosis, presence of ovarian follicles larger than 10 mm or other ovarian cystic lesions, chromosomal abnormalities, history of ovarian surgery, chemotherapy/radiotherapy, polycystic ovary syndrome, known autoimmune diseases, chronic, systemic, metabolic and endocrine pathologies, history of drug use. The sample size was estimated a priori using a two-sample test for proportions, according G*Power software to obtain >95% CI. All participants provided written informed consent. Genomic DNA was extracted from peripheral blood samples using the PureLink ® Genomic Kit (Invitrogen, Carlsbad, CA, USA) following the manufacturer’s protocol. DNA quality and quantity were assessed using the Qubit Fluorometer (Invitrogen) and the 4200 Tape Station Instrument (Agilent Technologies, Inc, Santa Clara, CA, USA). Only high-quality DNA samples were used for subsequent library preparation and sequencing. NGS was performed using the Ion AmpliSeq™ POI Panel on the Ion Torrent platform (Thermo Fisher Scientific, MA, USA). This custom targeted NGS panel provides complete exon coverage of 72 genes known to be associated with POI ( 34 ) (see details in Supplementary File S1 ). Library preparation was performed using10ng of genomic DNA measured with the Qubit 2.0 fluorometer (Thermo Fisher Scientific). Libraries were prepared both manually according to the Ion AmpliSeq Library Kit Plus and automatically on the Ion Chef™ instrument using the Ion AmpliSeq Kit for Chef DL8 (Thermofisher Scientific). Libraries were deemed suitable for sequencing after quality and quantity assessment. The concentration of each cDNA library, which was prepared manually, was determined on the Agilent 4200 system using the Agilent High Sensitivity DNA Assay (Agilent Technologies) according to the manufacturer’s recommendations. The concentration of the automatically prepared cDNA library pool was determined using the Ion Library TaqMan Quantitation Kit (Thermofisher Scientific). Libraries were diluted to 30 pM and then loaded into the Ion Chef™ instrument (Thermofisher Scientific). The Ion Chef™ instrument utilized the Ion 510™ & Ion 520™ & Ion 530™ Kit – Chef (Thermofisher Scientific) to perform emulsion PCR, enrichment and loading of the Ion S5–520 and/or 530 chip. PCR products obtained with the BigDye Direct Cycle Sequencing Kit (Applied Biosystems, Thermofisher Scientific) were purified with the BigDye Terminator Purification Kit (Applied Biosystems) and sequenced with the BigDye Terminator v3.1 (Applied Biosystems, Foster City, CA, USA) using the Genetic Analyzer 3500 Dx (Applied Biosystems). The exons of DNAH5 (ex 66), LAMC1 (ex 17), AR (ex 1), GDF9 (ex 3), ADAMTS19 (ex 2) were amplified using the primer list in Supplementary File S1 . The electropherograms were analyzed with SeqScape software, v.4 (Applied Biosystems) and compared with the corresponding reference genomes (DNAH5: NM_001369.3 ; LAMC1: NM_002293.4 ; AR: NM_000044.6 ; GDF9: NM_005260.7 ; ADAMTS19: NM_133638.6 ). The raw data generated from NGS were processed using Torrent Suite software v.5.14 which performed sequence alignment, adapter trimming, signal filtering, and quality-based read exclusion. Coverage analysis was performed using the Coverage Analysis plug-in. The coverage criteria for removing sequences from further analysis were: Reads < 200,000; mean depth < 200; uniformity < 90% (see Supplementary File S2 for details). Sequence variants detected within the 1549 amplicons of the 72 genes of the custom panel were analyzed with Ion Reporter version 5.18.2 using ‘Torrent Variant Caller v5.18-2 (Thermo Fisher Scientific). The parameters used to filter the variants were coverage, quality and frequency. The following parameters were set for the POI samples: Coverage ≥ 200, Quality score ≥ 30. To remove common germline variants, variants detected in the POI cohort were filtered against those found in the control group (N = 200). The POI candidate variants were further filtered by the non-Finnish European population from GnomAD ( https://gnomad.broadinstitute.org/ ) and the European population from 1000 Genomes ( https://www.internationalgenome.org/ ). The potential deleterious effects of the identified variants were assessed using the prediction algorithms SIFT and Polyphen2 ( 36 ), CADD ( 37 ), Grantham ( 38 ), MutatioTaster ( 39 ). Finally, the semi-automated process of the InterVar web service ( 40 ) and a manual review were used to identify the resulting potentially pathogenic variants. Both methods adhered to the American College of Medical Genetics and Genomics (ACMG) guidelines ( 41 ). The criteria used based on Richards et al. are listed in Supplementary File S3 . CNV analysis was performed using Ion Reporter software version 5.18.2. Briefly, a Hidden Markov Model (HMM)-based algorithm used normalized read coverage across amplicons to predict ploidy values (0, 1, 2, 3, etc.). Prior to CN determination, read coverage is corrected for GC bias and compared to a baseline coverage created from a control group of regions with known ploidy status. For each sample, a MAPD (Median Absolute Pairwise Difference) value is calculated, a metric that measures the noise in read coverage across all amplicons. To create a CNV call, MAPD is <0.4 was used as criteria. We considered the confidence score that filters out CNV regions that are likely to be false positives ( 42 ). As recommended by Ion Reporter software, we considered significant only CNVs that showed confidence score > 10. According to commonly used parameters ( 43 ), focal CNVs were defined as aberrations <3Mb in size. Three supervised learning-based tools (STRVCTVRE ( https://strvctvre.berkeley.edu/ ), CADD-SV ( https://cadd-sv.bihealth.org/ ) and X-CNV ( http://119.3.41.228/XCNV/search.php ) were used to predict the impact of CNV on pathogenicity. These predictors for the effect of structural variation (SV), CADD-SV ( 44 ) and StrVCTVRE ( 45 ), are trained with a random forest classifier and use a set of genomic features for SVs related to conservation, gene importance, coding region, expression and exon structure. X-CNV ( 46 ) is a framework based on the probabilistic value of the XGBoost algorithm that generates a meta-voting prediction (MVP) score to quantitatively measure the pathogenic effect of CNVs. As stated by the developers ( 44 ), the CADD-SV score on the Phred scale ranges from 0 (potentially benign) to 48 (potentially pathogenic). The StrVCTVRE score ranges from 0-1, with a score of 1 being more harmful ( 45 ). The MVP score of the X-CNV tool ( 46 ) indicates that CNVs between 0.46 and 0.76, between 0.16 and 0.46 and between 0.14 and 0.16 have potentially likely pathogenic, uncertain and likely benign effects, respectively. Genomic DNA was prepared using standard methods. Quantitative real-time PCR (Q-PCR) was performed using the Power SYBR Green PCR Master Mix with the QuantStudio 12K Flex Real Time System (ThermoFisher), as previously described ( 47 ). Normalization was performed to the GAPDH DNA content. Relative DNA amounts were calculated using the comparative cycle threshold method ( 48 ). Primer sequences are listed in Supplementary File S1 . Statistics was performed by Student’s t-test.