Targeted Next-Generation Sequencing Indicates a Frequent Oligogenic Involvement in Primary Ovarian Insufficiency Onset

Primary ovarian insufficiency (POI) is one of the major causes of female infertility associated with the premature loss of ovarian function in about 3.7% of women before the age of 40. This disorder is highly heterogeneous and can manifest with a wide range of clinical phenotypes, ranging from ovarian dysgenesis and primary amenorrhea to post-pubertal secondary amenorrhea, with elevated serum gonadotropins and hypoestrogenism. The ovarian defect still remains idiopathic in some cases; however, a strong genetic component has been demonstrated by the next-generation sequencing (NGS) approach of familiar and sporadic POI cases. As recent evidence suggested an oligogenic architecture for POI, we developed a target NGS panel with 295 genes including known candidates and novel genetic determinants potentially involved in POI pathogenesis. Sixty-four patients with early onset POI (range: 10–25 years) of our cohort have been screened with 90% of target coverage at 50×. Here, we report 48 analyzed patients with at least one genetic variant (75%) in the selected candidate genes. In particular, we found the following: 11/64 patients (17%) with two variants, 9/64 (14%) with three variants, 9/64 (14%) with four variants, 3/64 (5%) with five variants, and 2/64 (3%) with six variants. The most severe phenotypes were associated with either the major number of variations or a worse prediction in pathogenicity of variants. Bioinformatic gene ontology analysis identified the following major pathways likely affected by gene variants: 1) cell cycle, meiosis, and DNA repair; 2) extracellular matrix remodeling; 3) reproduction; 4) cell metabolism; 5) cell proliferation; 6) calcium homeostasis; 7) NOTCH signaling; 8) signal transduction; 9) WNT signaling; 10) cell death; and 11) ubiquitin modifications. Consistently, the identified pathways have been described in other studies dissecting the mechanisms of folliculogenesis in animal models of altered fertility. In conclusion, our results contribute to define POI as an oligogenic disease and suggest novel candidates to be investigated in patients with POI.


INTRODUCTION
Female factors account for one-third of all causes of infertility. Besides tubal disease and endometrial pathology, the dysregulation of any essential step involved in the ovulation of a competent oocyte may cause primary ovarian insufficiency (POI), a clinical syndrome defined by the premature loss of ovarian function. A recent meta-analysis of 31 epidemiological studies on the prevalence of POI in different countries between 1987 and 2018 reports an overall occurrence up to 3.7% in women younger than 40 years (1). At present, this disease is currently diagnosed when fertility is irreversibly affected. POI can manifest with a wide variety of clinical phenotypes, ranging from ovarian dysgenesis (OD) and primary amenorrhea (PA) to post-pubertal secondary amenorrhea (SA) for more than 4 months with raised gonadotrophins and low estradiol. The consequent long-standing estrogen deficiency exposes these women to an increased risk of complications such as cardiovascular diseases, reduced bone mineral density, and cognitive impairment (2). This disorder is highly heterogeneous in its etiology and several causes have been reported, mainly genetic, associated with chromosomal abnormalities (especially including X chromosome, such as in Turner syndrome), but also autoimmune, infectious, or iatrogenic. However, most causes of POI are still unknown, and the identification of novel causative genes is challenging. More recently, the advent of next-generation sequencing (NGS) technique and, especially, the whole exome screening (WES) of large POI families expanded the list of candidate genes to be screened in patients and consequently empowered the possibilities of a genetic diagnosis in idiopathic cases (3). Some WES studies demonstrated that pathogenic variants in meiotic chromosome pairing and synaptonemal complex (4)(5)(6) or alterations of other proteins of DNA recombination and repair (7,8) could be responsible for POI onset by usually impairing meiotic progression and triggering oocyte death, as further evidenced by murine models (9). Other studies identified variants in the folliculogenesis players of all stages of ovarian follicle maturation, which involves the precise interaction of hundreds of genes: from the primordial follicle stock establishment of ovarian reserve (10) to the primordial to primary follicle activation (11,12), throughout the follicular development in the gonadotropin-independent (13,14) and gonadotropin-dependent stages (15,16). Furthermore, alterations in genes contributing to extracellular matrix (ECM) remodeling by proteolytic activity on specific substrates within the ovarian context have been linked to the ECM turnover of abnormal somatic cells, thus leading to defects in follicular development (17). Since biological processes related to metabolism and immune system activation resulted to be enhanced in gene expression dynamics along with ovary development, including pathways associated with cell cycle, proliferation, apoptosis, ovulation, angiogenesis, and steroidogenesis (18), the disturbance of any of these pathways has been associated with reproductive diseases like POI (19)(20)(21)(22). Moreover, the remarkable point that emerged from recent NGS studies is the occurrence of oligogenic defects (6,23). From this perspective, various interacting genes might affect several mechanisms and pathways, and the synergistic and/or cumulative effect of several variants may contribute to POI phenotype. The NGS results in 64 patients with PA or early onset SA based on our panel of 295 POI candidate genes presented here aims to contribute to expand the list of potentially causative candidate genes and to support a genetically heterogeneous architecture of POI, resulting from defects in multiple complementary pathways.

Next-Generation Sequencing Panel Construction and Analysis
Genomic DNA (gDNA) was extracted from peripheral blood of enrolled patients. The Ampliseq Custom DNA panel (Illumina) was designed ad hoc to include the coding exons and flanking splice sites of 295 genes. The NGS panel, called OVO-Array (Table S1), has been obtained by joining data from: a) literature search of genes known to be previously associated with female infertility, with a role in ovary and in menopause onset (n = 159) (23); b) transcriptomic analysis of human granulosa cells treated with the oocyte-derived growth factor BMP15 (n of genes: 19; focus of another study with manuscript under review); and c) WES on 10 Italian women selected with the most severe phenotypes of PA and OD, six of them with familial inheritance (n = 117). The latter were identified among a total of 18,570 rare variants (missense, 86%; nonsense, 3.5%; indels, 5.6%). Briefly, pathway analysis (by Reactome v.74) of 1,916 genetic identifiers, found mutated at least once in WES, revealed an enrichment in the following biological processes with a key role in ovarian functionality (e.g., chromatin organization, cell cycle and meiosis, extracellular matrix organization, and cell-cell communication). Within these pathways, we further selected 117 genes potentially correlated with POI onset. The total coverage of the target genes by the designed amplicons was 100%. Library was prepared using enzymatic DNA fragmentation, with 50 ng of total gDNA and quantified with Quant-iT PicoGreen (Thermo Fisher Scientific). Nextera Rapid Capture Enrichment protocol (Illumina) was followed to tagment gDNA, amplify tagmented gDNA, hybridize probes, capture hybridized probes and for library capture and amplification. The library was then loaded onto the reagent cartridge (Illumina) and sequencing was performed on a NextSeq 500 (Illumina). Reads were then examined to identify single nucleotide alterations or small insertions/deletions, and a first in silico analysis (including base calling and demultiplexing) has been performed using MiSeq provided software (Real Time Analysis RTA v.1.18.54 and Casava v.1.8.2, Illumina). FastQ files for each sample, containing mate paired-end reads after demultiplexing and adapter removal, have been used as input for MiSeq pipeline. Briefly, FastQ files have been processed with MiSeq Reporter v2.0.26 using the Custom Amplicon workflow. This analytical method required FastQ files, a "Manifest file" containing information about the sequences of primer pairs, the expected sequence of the amplicons, and the coordinates of the reference genome (Homo sapiens, hg19, build 37.2) as input. Each read pair has been aligned using the MEM algorithm of the BWA software. Local Indel realignment and base recalibration step were performed using the software GATK. The realigned and recalibrated BAM file was used as input to GATK Unified Genotyper thus generating a VCFv4.2 file for each sample. Quality control of sequencing data was performed directly on FastQ files using the FastQC software. Reads were also filtered based on quality mapping and removed if their quality mapping was <20. Genetic variants showing a PHredScore lower than 20 were also filtered out.

POI Patients
We collected a cohort of patients to be analyzed through the OVO-Array panel, following IRCCS Istituto Auxologico Italiano Ethics Committee approval (BIOEFFECT, code 05M101_2014), composed of a total of 64 women with either primary (PA, n = 21) or early onset secondary amenorrhea (SA, n = 43), with onset before 25 years, all characterized by FSH values >40 IU/L and low estradiol. We excluded karyotype abnormalities, FMR1 premutations, and ovarian autoimmunity in all of them. Moreover, we selected additional 43 patients (PA, n = 18 and SA, n = 25) with the same inclusion criteria and early onset POI before 25 years of age. These patients have been analyzed between the years 2013 and 2020 for variants on a selected subset of only nine causative genes for POI, mainly involved in folliculogenesis and meiosis (BMP15, FIGLA, FOXL2, FSHR, GDF9, NOBOX, NR5A1, SYCE1, STAG3) for diagnostic purposes at IRCCS Istituto Auxologico Italiano. The results obtained with the different screening approaches (OVO-Array vs. diagnostic routine) have been then compared. Informed consent was obtained from all patients prior to blood sample collection and molecular studies.

In Silico Analysis of NGS Results and Variant Interpretation
The genetic variants resulting from the OVO-Array experiment and those derived from the diagnostic routine were annotated using Annovar software (24) and then classified as rare if resulting unknown or with a minor frequency allele (MAF) <0.01 in 1000 genomes, dbSNP, or EXAC databases. By checking for pathogenicity prediction the VarSome database until July 2021 (25), only those rare variants classified as likely pathogenic (LP), pathogenic (P), or variant of unknown significance (VUS), according to the American College of Medical Genetics (ACMG) classification guidelines (26), were considered for further analysis and confirmed by using Sanger sequencing. Gene ontology analysis of all the significatively altered genes was performed against DAVID Bioinformatics v6.8 (27) and Reactome Pathway Browser version 74 (28) on October 26, 2020. Both databases cross-reference other resources (e.g., NCBI, Ensembl, UniProt, KEGG, ChEBI, PubMed, and GO).

NGS Identified Novel Variants in Putative Candidates for POI
We recruited 64 patients presenting PA or SA with early onset (from 10 to 25 years) and performed NGS analysis on 295 selected genes, including known candidates and novel potentially causative genes. The mean coverage depth of the target regions was greater than 90% at 50× for all patients. Variants with MAF >0.01 were filtered out. We considered the identified variants relevant only if they were rare (MAF <0.01) or never described; predicted as P, LP, or VUS; and already known to be associated with POI. We also considered nine variants predicted as likely benign, found in recurring genes with a higher frequency in our patients than in the general female population or previously associated to POI and functionally characterized.
Moreover, among the 43 patients screened for diagnostic purposes, we could identify at least one genetic variant in known POI genes in 11 of them, thus providing a genetic diagnosis in about 25% of POI patients through NGS. Only one patient harbored two different alterations and another one harbored three variants ( Figures 1B, 2B, Tables S4, S5). All the identified variants were confirmed by Sanger sequencing.

Correlation Among the Number of Variants per Patient, Predicted Pathogenicity, and Phenotype
We summarized the prioritized variants harbored by each patient analyzed by OVO-Array NGS panel in Supplementary Tables S2 and S3. The phenotype of each patient is also reported. Nomenclature validation was performed using Mutalyzer 2.0.33 (https://www.mutalyzer.nl/) according to HGVS nomenclature version 2.0.
Tables 1-6 display all the genetic variations identified by OVO-Array by the number of variants per patient. Briefly, each table prioritizes patients depending on their phenotype severity and shows a detailed characterization (including origin, karyotype, age of menarche and POI onset, and familiarity). The frequency of each variant is then compared with those of the female general population available in the gnomAD ver.2.1.1 public dataset. The number of variants per patient is thus correlated with the VarSome predicted pathogenicity and phenotype of the patients. For each variation, VarSome criteria and its relative updated link are shown.

NGS Analysis Revealed Novel Rare Variants Affecting Pathways Involved in Ovarian Physiology
Taken together, our results uncovered 12 already known variants associated to POI: p.Arg300Leu in REC8; p.Arg68Trp in BMP15; p.Glu122Lysfs*45 in FIGLA; p.Pro103Ser, p.Thr121Ile, and p.Pro374Leu in GDF9; p.Phe543Serfs*7, p.Gly111Arg, and p.Lys371Thr in NOBOX; p.Val355Met in NR5A1; p.Arg643* in TP63; and p.Asn153His in LARS2 (green colored genetic variants in Tables S2, S4). Moreover, we could identify 41 novel rare variants in genes already associated to POI etiology (red colored in Tables S2, S4) and 74 novel rare variants in additional genes participating in key steps of ovarian follicle development, which may be considered putative candidates involved in POI pathogenesis (Table S2). Interestingly, 74 variants have been identified in 27 recurring genes among our two POI populations, and some of these genes/variants have already been associated to POI: 9 genetic alterations and 10 genes (green and red colored in Table 7, respectively). Furthermore, we found 12 alterations in recurrent genes which affected more than one patient and resulted at higher frequency in our POI populations, with respect to the general female population reported in the gnomAD (ver. 2.1.1) database ( Table 7)

Gene Ontology Analysis
The gene ontology analysis was carried out by interrogating the Reactome pathway database and DAVID gene functional annotation tool (refer to Tables S6, S7 for the statistical evaluation of the enrichment analysis). Considering both bioinformatic tools, we identified five main pathways likely affected by multiple genes found altered in our OVO-Array cohort and with functions in ovarian development and physiological maturation of oocytes ( Figure 1A). Most altered genes, 34 in total, participated in meiosis, cell cycle, and DNA repair processes, whereas 15 genes were involved in reproductive pathways (folliculogenesis, oocyte maturation, and follicular development). Both groups resulted to be affected in all manifestations of POI, from 46,XX OD to SA, but especially in the most severe phenotypes. Twelve genes with variations were identified in the ECM remodeling pathway, which influences relevant processes during follicle development, such as cell morphology, communication, proliferation, survival, and steroidogenesis, and appeared to be mostly affected with later POI onset (SA <25 years), but also in few exceptional 46,XX OD or PA cases. Fourteen altered genes were involved in the control of specific aspects of cell metabolism with roles in ovarian function, in both PA and SA patients. A rare frameshift alteration of SAMD11, which was demonstrated as a promoter of cell proliferation, has been found in several patients mostly affected with SA. Four different variants of RYR3, coding for a calcium channel with a role in the homeostasis of calcium, have been identified among all phenotypes. Furthermore, another identified pathway, fundamental in granulosa cell differentiation and proliferation, is the NOTCH signaling, represented by NOTCH2, NOTCH3, and NOTCH4 genes found to be altered in both PA and early SA cases. Finally, we also identified alterations in genes belonging to WNT signaling, cell death, and post-translational modifications mainly affecting SA patients (Table S2, Figure 3).
Most patients showed at least two pathways affected by variable genetic variations and a few showed several variations in genes from only one pathway (e.g., patient 27, presenting with PA and one variant in TP63 and FANCA genes, both related to meiosis, cell cycle, and DNA repair processes). Among SA, we also found three patients with the involvement of several variations in a singular pathway: patient 19 with three variants and patient 31 with two variants affecting genes related to meiosis, cell cycle, and DNA repair processes and patient 23 with a variant in BMP15 and a compound heterozygosity in GDF9, both key actors of folliculogenesis (Table S3).

DISCUSSION
In this study, we used the OVO-Array panel, a targeted NGS approach, to investigate the genetic cause of POI in a cohort of 64 patients with early onset of the disorder. The sequencing was focused on 295 candidate genes selected from literature and our previous data. Through this approach, we could identify 9 already known variants (6, 31, 34, 36-38, 42, 43) and 32 novel rare variants in genes already associated to POI etiology. Moreover, 34 novel rare variants have been found in additional genes participating in pathways important for ovarian physiology.      Our data are consistent with several lines of evidence recently emerging pointing to an oligogenic architecture for POI, according to which the presence of multiple genetic variants may partially explain the heterogenicity of phenotypes observed in POI patients. Indeed, we observed at least two variants in 35 patients, 5% of the cohort screened with our OVO-Array panel. Noteworthy, two of these patients harbored six variants. On the contrary, the screening of nine POI genes for diagnostic routine permitted to identify only 2 out of the 11 patients with alterations (18%) carrying either two or three variants. We observed that almost all patients screened by the OVO-Array, who carried the major number of variants, presented with a severe phenotype (PA or early SA onset just after menarche), whereas the majority of women carrying only one or two alterations (15 out of 25) were affected by SA. However, we also found that some patients harboring fewer pathogenic or likely pathogenic variants displayed a severe phenotype. These findings suggest that POI could emerge either by the disruption of a single fundamental genetic function or as a result of multilocus variations in genes interacting within different pathways, as proposed by other authors for oligogenic diseases (44). On this line, patient 6 of this cohort was previously reported for the potential synergic detrimental effect of a complex pattern of multiple inherited genetic variants in FIGLA, NOBOX, and NR5A1 (36). In the present study, we included this patient into the OVO-Array analysis, and by sequencing additional genes, we could identify another alteration (c.3755C>A, p.Pro1252Gln) in NCOR2. This gene encodes for a nuclear receptor corepressor that, together with its paralog NCOR1, mediates the retinoic acid-dependent repression of Fgf8 in mouse during organogenesis and, if this complex is mutated, exhibits increased Fgf8 expression, similar to retinoic acid deficiency (45). On one hand, retinoic acid is critical for the entrance in meiosis of ovarian germ cells (46); on the other hand, FGF8 cooperates with BMP15 to promote glycolysis in cumulus cells (47). Therefore, we could hypothesize that an alteration affecting these fundamental pathways might combine with the other three already described variants and play a synergic effect leading to the severe phenotype of the patient. This evidence supports multilocus analysis as a fundamental tool to explain the full phenotypic spectrum of women with POI.
Since NGS has demonstrated to be a powerful tool for the identification of new molecular players and pathways in POI onset, the information derived from the analysis of large NGS panels, such as the OVO-Array, might increase the diagnostic power up to 75% of POI cases, in contrast to the current 25% of positive diagnosis obtained by screening few POI genes. The actual limit of this approach is the necessity of validating new putative candidates by further functional testing and sequencing in larger cohorts. Nonetheless, in this study, we could observe both the presence of the same variant in more than one patient at an increased frequency in women with POI with respect to the general female population and different alterations affecting the same genes in several patients. Although the fertility history of      the controls reported in public datasets is unknown, thus introducing a potential bias in the analysis, our findings support the relevance of these candidates in POI pathogenesis and deserve further investigations.

Large-Scale Genetic Analysis Identified Novel Variants in Genes Participating in Pathways Relevant for Ovarian Function
A first hint pointing at the possible involvement of the newly identified candidates derives from gene ontology, which highlighted pathways fundamental for the reproductive success in mammals. Meiotic recombination is a complex process, which requires the combination of a large plethora of factors (48). Furthermore, during the arrest at meiotic prophase I, primordial follicle oocytes are more vulnerable to DNA double stand breaks emanating from endogenous and exogenous sources (49); thus, pathogenic variations in one or more of their encoding genes might explain the onset of POI, as previously proposed (50). Noteworthy, we found 34 altered genes participating in meiosis, cell cycle, and DNA repair as possible candidates for POI. This pathway resulted to be mainly affected in patients with the most severe phenotypes, and in several cases, we found variants in at least two genes encoding DNA repair and meiotic factors. Moreover, this was the only pathway affected by genetic variants in some patients. Our data indicate that, in line with its oligogenic nature, POI might also result from alterations in more loci belonging to the same pathway.
The next major pathway in which we have found several altered genes was folliculogenesis, encompassing all stages of ovarian follicle development. In one case of SA (patient 23), we found variations in two major factors of folliculogenesis: a variant of BMP15 (p.Arg68Trp) with deleterious functional effects (31), together with a compound heterozygosity in GDF9 [p.Tyr93Cys and the previously reported p.Pro374Leu (34)], supporting the possibility that the disruption of a unique pathway at variable levels can also cause the disorder. Besides the known actors (i.e., FIGLA, NOBOX, BMP15, GDF9, FSHR, LHCGR, NR5A1, AR), we found variants in two genes involved in autophagosome assembly (ATG2A and ATG4C). Their role may deserve further attention because autophagy is an important mechanism in mammalian ovarian development, regulating follicular atresia (51). We also identified genes causing disorders of sex development or with roles in spermatogenesis. Indeed, CYP21A2 alterations are responsible for congenital adrenal hyperplasia, caused by diminished aldosterone and cortisol production that result in ambiguous genitalia in female affected (52). DMRT3 encodes for a transcription factor with evolutionary conserved roles in sex development, whose haploinsufficiency leads to 46,XY male to female sex reversal (53). TUBA8 encodes an isoform of a-tubulin highly expressed in ovarian follicle (www.proteinatlas.org). This gene has an ortholog in mouse which is expressed in the brain and testis, with a role in spermatid development (54), and can cause asthenozoospermia in men (55). ID1 has been described as potential AMH downstream target genes (56,57), and since it is involved in the regulation of follicular growth, further characterization of its molecular function in this pathway would be needed. KMT2D alterations are the main cause of Kabuki syndrome, a congenital intellectual disability with genitourinary anomalies among other additional features. Interestingly, a reduced number of dominant families was reported in a Kabuki cohort (58), thus supporting a possible role of KMT2D variations in female fertility.
We also identified alterations in genes involved in ECM remodeling in the OVO-Array POI cohort, such as ADAMTS4, ADAMTS5, ADAMTS16, AGRN, COL6A1, COL6A2, ERBB3, ERBB4, PKP1, RELN, THBS2, and VWF. Given the ECM FIGURE 3 | Bubble chart of the pathways identified by the gene ontology analysis. The sample numerosity (n) is reported for each phenotype (46,XX OD, light red; PA, light blue; eSA, light yellow; and SA <25 years, light green). Pathways are represented by colored bubbles (refer to the legend on the right). The size of the bubbles depends on the number of variants identified in each pathway for each phenotype. The percentages within the bubbles indicate the ratio between the number of variants per pathway associated to a phenotype and the total of variants associated to that phenotype. 46,XX OD = ovarian dysgenesis with normal karyotype (46,XX); PA = primary amenorrhea with functional ovarian defect, eSA = secondary amenorrhea onset few months after menarche; SA = secondary amenorrhea <25 years of age.
contribution to granulosa cell survival and proliferation, the ECM is required for maintaining the follicular cell morphology in all phases of ovarian follicle function (59). This ECM role is supposed to be mediated, at least in part, by the distinct ADAMTS subtypes and collagens. The ECM regulates cell aggregation and intracellular communication between the oocyte, granulosa, and theca cells within the follicle. The communication of small metabolites, ions, and second messengers between each cell type is made possible by a network of gap junctions, such as desmosomes containing among others the accessory plaque protein plakophilin 1. The ECM proteins are also necessary to support granulosa cell survival and steroidogenesis; for example, reelin contributes to follicular stability. Factors within the ECM are important in the control of follicular development and atresia. Included in this group of ECM proteins are thrombospondins, which are thought to contribute to the regulation of angiogenesis, and von Willebrand factor (59).
Several elements of Notch signaling appear to be involved in the POI pathogenesis in the OVO-Array cohort (NOTCH2, NOTCH3, NOTCH4). The Notch pathway is a contactdependent signaling system active in the mammalian developing ovary which has multiple functions in follicle assembly, maturation, development, and meiotic entry (60). Notch proteins (NOTCH1, NOTCH2, NOTCH3, and NOTCH4) function as transmembrane receptors for one of the membrane-bound ligands (JAG1, JAG2, DLL1, DLL3, and DLL4) and their binding causes a conformational change of the Notch protein starting a series of sequential proteolytic cleavages at the receptor juxtamembrane region allowing the release of the Notch intracellular domain that is eventually free to translocate to the nucleus. Within the nucleus, it activates the transcription of Notch target genes. Receptors, ligands, modulators, and activated genes belonging to this system are expressed and finely regulated during folliculogenesis. Although Notch1 and Notch4 expression is limited to the ovarian vasculature, Notch2 and Notch3 are expressed in the granulosa cells of developing follicles and mediate follicle assembly and growth in mammals. Consistent with its role within the mammalian ovarian follicle, genetic alterations in NOTCH2 were previously identified via whole-exome sequencing and functional evidence demonstrated the correlation of NOTCH2 missense variations in POI (23,61).
Several rare variants have been identified in genes regulating different aspects of follicular cell metabolism. Hexokinase 3 (HK3) converts glucose to glucose-6-phosphate in the first step of glucose metabolism, exerting a protective effect against oxidative stress (62), and it was associated with age at natural menopause (63). A recent transcriptomic study in primates demonstrated that genes most regulated in aged oocytes and granulosa cells are related to antioxidant defenses (64), further supporting the importance of cellular damage in infertility. POLG and LARS2 play roles in mitochondrial DNA replication, gene expression, and protein synthesis and degradation (65). Mutations in POLG can cause among others the neurological conditions Alpers syndrome and progressive external ophthalmoplegia (PEO), and the latter can associate with POI (35). Four POLG variants emerged by our analysis: a PA (patient 23, see 3.1.1 section) and a SA patient with c. c.752C>T (p.Thr251Ile) and c.1760C>T (p.Pro587Leu) in compound heterozygosity; c.803G>C (p.Gly268Ala) and c.3436C>T (p.Arg1146Cys) in two patients with SA. Unfortunately, no further neurological information was available for these patients. Mutations in LARS2, encoding mitochondrial leucyl-tRNA synthetase, lead to POI and hearing loss in Perrault syndrome (66). We identified the c.457A>C (p.Asn153His) variant (ClinVar VCV000191173.2) of LARS2 in a patient with SA, who was not previously associated with Perrault syndrome, and other three missense variants in two SA patients. Unfortunately, we have no further clinical information also in these patients. The cholesterol synthetase DHCR24 belongs to the pathway of steroid biosynthesis which regulates many physiological processes (i.e., stress response, ovarian cycle, and endocrine system) (67), and this is the first report of a variant in DHCR24 in a patient with POI. The VLDLR gene in humans is relevant for steroidogenesis as its expression in granulosa cells of pre-ovulatory follicles keeps a relevant role in lipoprotein endocytosis during follicular growth (19). Evidence in hens correlates the VLDLR function with fertility (68) and its haploinsufficiency in women has been associated to impaired folliculogenesis (69). Here, we confirm previous findings by describing a novel missense variant in a patient with PA.
WNT signaling has essential functions in ovary differentiation, follicle development, and hormone synthesis (70). In this study, we identified missense variants in LGR4, one of the receptors for R-spondins, which augment the WNT signaling pathway (71), and LRP5, a co-receptor for the canonical Wnt-b-catenin signaling with a known role in the regulation of bone mineral density (72), but we could not find any further evidence in the ovary.
Cell proliferation and death pathways were represented by the identification in our cohort of variants in SAMD11 and RIPK1, respectively. SAMD11 has never been associated with ovarian phenotypes before. This gene encodes for a transcriptional modulator relatively uncharacterized but phylogenetically conserved from zebrafish to humans, widely expressed in various tissues (73), where it might promote cell proliferation. Conversely, RIPK1 can mediate apoptosis during embryonic development (74). This is the first report of variants in these genes associated with an ovarian phenotype, and their role in the ovary should be investigated in further studies. Another gene involved in mitosis control that we found altered in POI was USP35, encoding a deubiquitinating enzyme which regulates the stability and function of Aurora B kinase and whose depletion results in inhibited metaphase chromosome alignment (75).
In mammalian oocytes, calcium is one of the major signal molecules involved in meiotic cell cycle resumption, arrest, and apoptosis. RYR3 encodes a calcium channel that mediates the intracellular release of Ca 2+ , and the ryanodine receptors (RyR) family has been identified in mammalian oocytes (76). Here, we reported several RYR3 variants with predicted pathogenic effect, which might dysregulate the homeostasis of calcium within oocytes, thus supporting a role for this gene in the disorder.

CONCLUSIONS
We propose novel candidate gene-disease variants likely causative or conferring susceptibility to POI onset. Our findings, by supporting the oligogenic nature of POI, suggest that in the presence of the most severe forms of ovarian insufficiency, it is necessary to screen multiple candidates or to perform WES analysis, since more alterations in different genes may synergize for the determination of the phenotype. Moreover, we show that multilocus analysis could increase the diagnostic power and the accuracy of POI diagnosis, thus ameliorating genetic counseling in patients. The systematic application of the OVO-Array multilocus analysis shall improve the management of POI including a personalized approach to the fertility defect and to the associated extra-ovarian abnormalities that can frequently anticipate or follow the POI onset.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material and online at https://doi. org/10.5281/zenodo.4543275 (Digital Object Identifier 10.5281/zenodo.4543275). Further inquiries can be directed to the corresponding authors.

ETHICS STATEMENT
The studies involving human participants were reviewed and approved by IRCCS Istituto Auxologico Italiano Ethics Committee. Written informed consent to participate in this study was provided by the patients or, if minors, by the legal guardian/next of kin of the participants.

AUTHOR CONTRIBUTIONS
RR was responsible for the conception and design, OVO-Array experiment, and interpretation of data and drafted the manuscript. SM was responsible for the sequencing analysis, carried out the in silico predictions, and helped in drafting the manuscript. FG was responsible for NGS for diagnostic screening. DG was responsible for the bioinformatic NGS analysis. LL contributed to the OVO-Array experiment. AM, CM, FB, and MB were responsible for the enrollment of patients and critical revision of the manuscript. LP designed and supervised the research and revised the manuscript critically for important intellectual content. All authors contributed to the article and approved the submitted version.