Oligogenic Inheritance Underlying Incomplete Penetrance of PROKR2 Mutations in Hypogonadotropic Hypogonadism

The role of the prokineticin 2 pathway in human reproduction, olfactory bulb morphogenesis, and gonadotropin-releasing hormone secretion is well established. Recent studies have highlighted the implication of di/oligogenic inheritance in this disorder. In the present study, we aimed to identify the genetic mechanisms that could explain incomplete penetrance in hypogonadotropic hypogonadism (HH). This study involved two unrelated Tunisian patients with HH, which was triggered by identifying a homozygous p.(Pro290Ser) mutation in the PROKR2 gene in a girl (HH1) with Kallmann syndrome (KS). The functional effect of this variant has previously been well demonstrated. Unexpectedly, her unaffected father (HH1P) and brother (HH1F) also carried this genetic variation at a homozygous state. In the second family, we identified a heterozygous p.(Lys205del) mutation in PROKR2, both in a male patient with normosmic idiopathic IHH (HH12) and his asymptomatic mother. Whole-exome sequencing in the three HH1 family members allowed the identification of additional variants in the prioritized genes. We then carried out digenic combination predictions using the oligogenic resource for variant analysis (ORVAL) software. For HH1, we found the highest number of disease-causing variant pairs. Notably, a CCDC141 variant (c.2803C > T) was involved in 18 pathogenic digenic combinations. The CCDC141 variant acts in an autosomal recessive inheritance mode, based on the digenic effect prediction data. For the second patient (HH12), prediction by ORVAL allowed the identification of an interesting pathogenic digenic combination between DUSP6 and SEMA7A genes, predicted as “dual molecular diagnosis.” The SEMA7A variant p.(Glu436Lys) is novel and predicted as a VUS by Varsome. Sanger validation revealed the absence of this variant in the healthy mother. We hypothesize that disease expression in HH12 could be induced by the digenic transmission of the SEMA7A and DUSP6 variants or a monogenic inheritance involving only the SEMA7A VUS if further functional assays allow its reclassification into pathogenic. Our findings confirm that homozygous loss-of-function genetic variations are insufficient to cause KS, and that oligogenism is most likely the main transmission mode involved in Congenital Hypogonadotropic Hypogonadism.

The role of the prokineticin 2 pathway in human reproduction, olfactory bulb morphogenesis, and gonadotropin-releasing hormone secretion is well established. Recent studies have highlighted the implication of di/oligogenic inheritance in this disorder. In the present study, we aimed to identify the genetic mechanisms that could explain incomplete penetrance in hypogonadotropic hypogonadism (HH). This study involved two unrelated Tunisian patients with HH, which was triggered by identifying a homozygous p.(Pro290Ser) mutation in the PROKR2 gene in a girl (HH1) with Kallmann syndrome (KS). The functional effect of this variant has previously been well demonstrated. Unexpectedly, her unaffected father (HH1P) and brother (HH1F) also carried this genetic variation at a homozygous state. In the second family, we identified a heterozygous p.(Lys205del) mutation in PROKR2, both in a male patient with normosmic idiopathic IHH (HH12) and his asymptomatic mother. Whole-exome sequencing in the three HH1 family members allowed the identification of additional variants in the prioritized genes. We then carried out digenic combination predictions using the oligogenic resource for variant analysis (ORVAL) software. For HH1, we found the highest number of disease-causing variant pairs. Notably, a CCDC141 variant (c.2803C > T) was involved in 18 pathogenic digenic combinations. The CCDC141 variant acts in an autosomal recessive inheritance mode, based on the digenic effect prediction data. For the second patient (HH12), prediction by ORVAL allowed the identification of an interesting pathogenic digenic combination between DUSP6 and SEMA7A genes, predicted as "dual molecular diagnosis." The SEMA7A variant p.(Glu436Lys) is novel and predicted as a VUS by Varsome. Sanger validation revealed the absence of this variant in the healthy mother. We hypothesize that disease expression in HH12 could be induced by the digenic transmission of the SEMA7A and DUSP6 variants or a monogenic inheritance involving only the SEMA7A

INTRODUCTION
Idiopathic hypogonadotropic hypogonadism (IHH) (MIM ID #146110) is defined as the absence of sex steroid synthesis associated with the lack of appropriate gonadotropin secretion. This leads to a variable degree of impuberism, often diagnosed during childhood or adolescence. When the embryonic migration of gonadotropin-releasing hormone (GnRH) neurons from the nasal placode to their final destination in the hypothalamus is disrupted, the resulting phenotype is Kallmann syndrome (KS), which is defined as the association of HH with hyposmia or anosmia (Sonne and Lopez-Ojeda, 2021). Nineteen genes are known to be involved in KS (ANOS1, FGF8, FGFR1, FGF17, CHD7, IL17RD, DUSP6, SPRY4, FLRT3, KLB, SEMA3A, SEMA3E, NSMF, HS6ST1, WDR11, SOX10, FEZF1, IGSF10, PROK2, and PROKR2) (Boehm et al., 2015;Stamou and Georgopoulos, 2018). Normosmic idiopathic hypogonadotropic hypogonadism (nIHH), which is not associated with anosmia, and results from the dysfunction of the GnRH neurons that successfully completed their embryonic migration to the hypothalamus. As of today, 46 genes have been associated with IHH (Topaloglu, 2017;Cangiano et al., 2020). The molecular alterations identified in these genes account for 40% of all IHH cases (Gach et al., 2020), thus suggesting that half of the IHH causal genes remain unknown. Molecular alterations have been identified for several neuropeptides or their corresponding receptors, which are involved in the physiology of the gonadotropic axis: GNRHR/GNRH1, KISS1R/GPR54, TAC3/TACR3, and PROK2/PROKR2 (Brioude et al., 2010;Topaloglu and Kotan, 2010). In addition to reproductive dysfunction, nIHH/KS patients may also manifest a variety of other non-reproductive disorders, such as midline facial defects, dental agenesis, renal agenesis, hearing loss, or bimanual synkinesis (Young et al., 2019). IHH may be inherited in an X-linked recessive, autosomal dominant, or autosomal recessive modes of inheritance. IHH has been predominantly detected in sporadic cases (Quaynor et al., 2011;Gach et al., 2020).
The PROK2 gene (NG_008275.1) is located on chromosome 3p21.1 and comprises four exons. The PROKR2 gene (NG_008132.1) maps to chromosome 20p13 and contains seven transmembrane-domain receptors. Both genes, PROK2 and PROKR2, belong to the family of prokineticins, and a group of multifunctional secreted proteins first identified in 2000 by Li et al. (2001). Their involvement in KS was strongly suggested in 2006 when homozygous Prokr2 knockout mice were shown to have hypogonadotropic hypogonadism due to hypothalamic GnRH deficiency and agenesis or hypoplasia of the olfactory bulbs (Matsumoto et al., 2006). This was confirmed in the same year by the discovery of genetic variations in patients with KS (Dodé et al., 2006). More than 27 genetic variations in PROKR2 (Dode and Rondard, 2013) and over 11 genetic variations in PROK2 have been reported in patients with nIHH or KS (Martin et al., 2011). However, the KS and nIHH genetics are complex and still not well understood (Dode and Hardelin, 2010). In addition, genetic variations in PROK2 and PROKR2 genes reported in KS and nIHH patients were found at heterozygous, homozygous, and compound heterozygous states (Dodé et al., 2006;Hardelin and Dode, 2008). Homozygous loss-of-function genetic variations in the PROK2 gene cause nIHH in mice and humans (Abreu et al., 2008;Schöneberg and Liebscher, 2021). Thus, an autosomal recessive mode of transmission was established (Pitteloud et al., 2007). Later, other studies on different ethnic populations reported a large number of heterozygous genetic variations in PROK2 and PROKR2 genes with considerable clinical and molecular heterogeneity among several patients having both KS and nIHH (Cole et al., 2008;Sarfati et al., 2010a). Heterozygous genetic variations inherited from clinically unaffected first-degree relatives in the PROKR2 gene have been reported in some KS/nIHH patients, which was attributed to digenic/oligogenic transmission rather than a dominant negative effect of monoallelic PROKR2 genetic variations (Dodé et al., 2006;Abreu et al., 2008;Cole et al., 2008;Monnier et al., 2009;Quaynor et al., 2011;Lewkowitz-Shpuntoff et al., 2012). The hypothesis of oligogenic inheritance postulates that disease expression is induced by more than one mutated IHH gene (Pitteloud et al., 2005;Maione et al., 2018). Indeed the PROKR2 gene has been involved in several digenic and trigenic associations such as PROK2/PROKR2, FGFR1/PROKR2, PROK/GNRHR, and PROKR2/CHD7/FEZF1 (Cole et al., 2008;Canto et al., 2009;Sarfati et al., 2010b;Pablo Méndez et al., 2015;Zhang et al., 2020). However, the expression of deleterious alleles is considerably variable if we compare the phenotypes of patients carrying identical variations (Pitteloud et al., 2007;Boehm et al., 2015). Several studies conducted on large cohorts have shown that oligogenic heredity accounts for 2.5-15% of Congenital Hypogonadotropic Hypogonadism (CHH) patients (Cassatella et al., 2018).
Here, we report a clinical and genetic investigation of KS and nIHH in two Tunisian families after excluding the involvement of monogenic inheritance of PROKR2 gene variants.

Patients
The current molecular analysis was conducted for two unrelated HH families (HH1-a woman with KS and HH12-a man with nIHH). The index cases were recruited at the Endocrinology Department of the National Institute of Nutrition in Tunis. Both patients belong to simplex families (only one family member is affected, with healthy siblings). For the two patients, GnRH deficiency diagnosis was established based on puberty state (absent or incomplete), hormonal tests (testosterone in HH12 and serum gonadotropins levels), response to GnRH, and anterior pituitary function, which was evaluated by measuring the basal levels of free T4, TSH, and prolactin as well as the cortisol peak levels after ACTH stimulation. Pituitary imaging was performed by magnetic resonance imaging (MRI) to exclude acquired causes of nIHH. Olfactory testing for anosmia or hyposmia was also assessed.

Molecular Investigation
Genomic DNA was extracted from peripheral blood leucocytes by FlexiGene DNA extraction kit (Qiagen) according to the manufacturer's instructions. The coding region and intron-exon boundaries of PROK2 and PROKR2 genes as well as exon 13 of the SEMA7A gene were amplified as previously reported (Dodé et al., 2006). The PCR products were sequenced with the Big Dye terminator kit (Applied Biosystems, Foster City, CA, United States) using one of the PCR primers on an ABI prism 3,130 DNA Genetic Analyzer (Applied Biosystems) in accordance with the recommendations of the manufacturer.

In-silico Analysis
To gain insight into the effect of the newly identified genetic variation p.(Lys205del), we carried out an in silico analysis. We first tested the impact of genetic variation on splicing by creating or abolishing a splice site. This analysis was performed using the Human Splicing Finder program (Desmet et al., 2009) 1 . We then carried out an in silico 3D structure prediction using Phre2 tool (Kelley et al., 2009) 2 in order to predict the effect of genetic variation on protein folding.

Exome Analysis
An in-house pipeline analysis was used to generate VCF files. Then, the annotation and the prioritization of variants were carried out using VarAFT tool, version 2.04 3 . A disease-causing gene approach was used to extract variants located in 46 CHHrelated genes extracted from the Online Mendelian Inheritance in Men database and recent literature (Supplementary Table 1). Functionally relevant variants located in exonic genomic regions and splice sites were then selected from the list of variants contained in the VCF file. Variants located in the regulatory regions (UTR, promoters, enhancers, and miRNA binding sites, etc.) flanking the PROKR2 gene were also screened. Regulatory regions were retrieved using the Genome Browser tool provided by the University of California Santa Cruz database 4 as well as the GeneCards Human Gene database 5 . The functional effects of the genetic sequence variants were evaluated by in silico prediction tools including SIFT 6 , PolyPhen 7 , FATHMM 8 , MutationTaster 9 , MutationAssessor 1.0 10 , PROVEAN v1.1 11 , and Varsome 12 . Furthermore, in order to pinpoint the variants that could potentially be associated with disease expression in the index cases (HH1 and HH12), each selected variant of interest was compared to the exome data of 92 unrelated Tunisian individuals stored in the local database. This approach allowed us to avoid bias especially that public databases, such as ESP6500, GnomAD, or 1000Genomes, do not contain relevant information regarding the genetic background of individuals from the Middle East, North Africa populations.

Digenic Effect Prediction
The Oligogenic Resource for Variant Analysis (ORVAL) online software 13 was used to predict the effect of variant combinations in disease-causing genes. This tool is based on variant annotation and effect prediction of two predictive methods named VarCOPP and Digenic Effect predictor. The evaluation scores of digenic combinations include the support score (SS) and the classification score (CS). SS informs about the percentage of algorithms integrated in VarCOPP and Digenic Effect predictor that support the pathogenicity of a given bi-locus combination. Its value ranges from 0 to 100. For candidate disease-causing combinations, SS should be equal or greater than 50%. The CS should correspond to the median probability of a digenic combination to be disease-causing. Its value ranges from 0 to 1. Variant pairs are considered pathogenic when CS is greater than 0.489. To estimate the probability that a predicted diseasecausing combination is a true positive, VarCOPP encompasses 95 and 99% confidence intervals which are delimited by minimal CS and SS values. The 95% confidence interval requires CS and SS values greater than or equal to 0.55 and 75, respectively. A digenic combination falls inside the 99% confidence interval when CS is greater than 0.74 and SS is equal to 100.
The variant combinations were categorized into three classes. The first one, termed "true digenic" involves two variants in two different genes to induce disease expression. Otherwise, if an individual carries only one variant, he will be considered unaffected. The second class is referred as "composite class." In this case, an individual carrying only one variant expresses the disease, whereas the second genetic locus is a genetic modifier that modifies the severity of the clinical presentation or the age at onset. The third class, called "dual molecular diagnosis, " requires the presence of variants in two different genes inducing two independent clinical entities (Papadimitriou et al., 2019).
Only the pathogenic bi-locus combinations for each individual were extracted from the list of variant pairs. Considering that variants found in common between the three HH1 family members would not be at the origin of disease expression, the digenic combinations that were specific to the index case were selected. The specificity criteria were related to the presence or absence of a variant or to its zygosity state, i.e., variants that are homozygous in the index case and heterozygous in the two asymptomatic family members.

Clinical Investigation
HH1 is a 21-year-old woman born to first-cousin phenotypically normal parents (Figure 1) and is originating from Central Tunisia. The patient presented congenital anosmia associated with obesity and absence of spontaneous puberty. Her basal and stimulated GnRH-gonadotropin levels were flat, and her MRI of the hypothalamic-pituitary region showed an aplasia of the olfactory bulbs ( Table 1). The 17-beta estradiol levels were not measured for this patient. Her family presents a history of lung cancer, congenital deafness, and colon polyps. Except for the cousin of the index case who had delayed puberty, the other family members had normal pubertal development and olfactory tests. HH12 is a 28-yearold man born to unrelated phenotypically normal parents (Figure 1). The patient suffers from nIHH associated with hypertension and obesity. He first consulted for impuberism. A clinical examination showed that his height is 1.78 m and his weight is 113 kg. The genitalia examination showed a micropenis with a gonadal volume of 3 ml. His GnRH response showed a normal response of follicle-stimulating hormone and luteinizing hormone. The testosterone serum level was 1.9 nmol/l (10-41.5 nmol/l). He presents a normal sense of smell and a normal MRI.  For the HH1 family, after having identified a missense genetic variation, c.868C > T; p.(Pro290Ser), in the PROKR2 gene at a homozygous state in the proband, we sequenced this gene in the relatives: the mother and the sister carried the p.(Pro290Ser) genetic variation at a heterozygous state, but, intriguingly, the father and the brother were homozygous. To confirm these molecular findings, blood sampling, PCR, and sequencing reactions were performed twice. To exclude contamination and paternity problems, all members of this family were genotyped using the identifier kit (Applied Biosystems). The analysis of allele segregation was in favor of paternity inclusion. For HH12, the molecular analysis allowed the identification of a novel variation, c.614_616del; p.(Lys205del), in the PROKR2 gene at a heterozygous state. Bioinformatics analysis using the Uniprot database showed that the lysine amino acid at position 205 is buried amid the second extracellular loop and that this residue is highly conserved in Prokr2 from mouse, rat, chimp, dog, cow, Xenopus tropicalis, and zebrafish and in the human PROKR1 protein. The deleterious effect of this genetic variation could be explained by the creation of a new acceptor splice site as illustrated by the results of the Human Splicing Finder and Genetic Variation Taster tools (donor gained, score = 0.83mut/=0.55 C). The prediction of the p.(Lys205del) variation effect performed with Phyre2 showed a deleterious effect on the 3D protein structure. The comparison between normal and mutant 3D structures was consistent with the deleterious effect. Indeed the lysine amino acid at position 205 interacts with other residues in the EC2, EC3, the third transmembrane domain (TM3), and the TM; thus, its deletion could reduce the PROKR2 protein stability (Figure 2). The genetic investigation of the mother of HH12 showed that she was a carrier of the same likely pathogenic mutation [p.(Lys205del)] at a heterozygous state.

Exome Analysis and Digenic Effect Prediction for the HH1 Family
In total, 43, 39, and 33 genetic variants were identified in HH1, her brother, and her father, respectively. For HH1, seven rare variants were identified with a minor allele frequency (MAF) < 0.1 according to the 1000Genomes and GnomAD databases. Four rare variants were harbored by the brother and the father. Besides the p.(Pro290Ser) in the PROKR2 gene, the remaining rare variants in each case were predicted as benign or as VUS by Varsome search engine as well as 10 in silico prediction programs (Supplementary Tables 2-4). For each set of variants identified in the three family members, we proceeded to the digenic effect prediction of variant pairs using the ORVAL software. A total of 940, 828, and 668 digenic combinations were, respectively, obtained in the index case HH1, the asymptomatic brother HH1F, and the father HH1P. In the proband, 62 variant pairs were predicted as pathogenic (Supplementary  Table 7) disease-causing combinations were found in the brother and father of HH1, respectively (Figure 3). The confidence intervals supporting the pathogenicity of the digenic combinations ranged from 90 to 99%.
The three family members carried 19 variants in 10 prioritized genes: PROKR2, DUSP6, NTN1, ANOS1, KISS1R, PROP1, DCC, IGSF10, SRA1, and PLXNA1. They shared almost 30 pathogenic digenic combinations involving these variants. In the case of HH1, the PROKR2 variant [p.(Pro290Ser)] yielded the highest number of bi-locus combinations (31%) predicted as diseasecausing, with 25 genes and a median pathogenicity score  of 0.62. We also found 13 pathogenic digenic combinations (Figure 3) involving variants that were either only carried by HH1 or those that were homozygous in HH1 and heterozygous in the asymptomatic cases ( Table 2). These specific variants included four exonic missense variants in PCSK1 (NM_000439), IL17RD (NM_017563), and CCDC141 (NM_173648) genes, which were homozygous in the proband and heterozygous in the unaffected family members. Another exonic variant in the FLRT3 gene (NM_013281), as well as two splice site variants in WDR11 (NM_018117) and SMCHD1 (NM_015295) genes, was only present in the index case HH1 (  (Figure 4). The variant pair CCDC141 (c.2803C > T)-PROKR2 (c.868C > T) was classified by ORVAL as true digenic. The contribution of the three missense variants in IL17RD and PCSK1 genes, which were homozygous in the index case and heterozygous in the asymptomatic cases ( Table 2), to the total number of pathogenic digenic combinations did not differ among the three family members (Figure 4). This indicates their minor contribution to disease expression in HH1. The splice site variant (c.5476 + 10A > G) in the SMCDH1 gene (NM_015295), only present in the index case ( Table 2), was not involved in any disease-causing digenic combination. The two other genetic variations identified in FLRT3 (NM_013281) and WDR11 (NM_018117) genes, present only in HH1 case, represented only 2% of the overall diseasecausing variant pairs.
We also evaluated the involvement of other variants located in the flanking regulatory regions of the PROKR2 gene, which could be responsible for the incomplete penetrance in the family members of HH1. On the basis of predictions provided by the Encyclopedia of DNA Elements and the GeneHancer database, we focused on regulatory regions adjacent to the PROKR2 gene, likely covered by whole-exome sequencing (WES). The closest promoter and enhancer regions are located in the exon 1-intron 1 junction of PROKR2. This region is rich with CpG islands and overlaps with the AX746654 gene. Two transcription factors, namely, HDAC2 and CTCF, are known to bind to this region. The WES analysis did not reveal any variants in the regulatory region of PROKR2 nor in the genes encoding its transcription factors.   (Glu587Lys)] was only present in HH12 and absent in his asymptomatic mother (Figure 1). The variants located in the promoter region of PROKR2 were extracted, which revealed one common variant (c.−9 + 342A > G) in intron 1 with a MAF of 0.3 according to GnomAD. Varsome is the only prediction tool that annotated this variant as benign. We have also identified in the CTCF gene, predicted as a transcription factor of the PROKR2 gene, one variant (c. * 29T > G) which is also frequent (MAF = 0.1) and predicted as benign by Varsome. Furthermore, the ORVAL prediction tool was employed to evaluate the digenic effect of both variants in combination with the p.(Lys205del) identified in PROKR2 gene. However, no digenic combinations were identified.

DISCUSSION
Since 2006, many studies have reported PROK2/PROKR2 genetic variations mainly missense, nonsense, and frameshift genetic alterations (Figure 5; Sarfati et al., 2010b;Sbai et al., 2014). Most patients carry monoallelic genetic variations, but some patients may carry bi-allelic variants either in the same or in different genes, indicating a dominant, recessive, or oligogenic inheritance, respectively (Dodé et al., 2006;Pitteloud et al., 2007;Leroy et al., 2008;Sarfati et al., 2010b). The p.(Pro290Ser) variant in PROKR2 gene has been described in CHH patients at heterozygous (Dodé et al., 2006) and homozygous (Sarfati et al., 2010a) states and is always associated with clinical presentations similar to those seen in HH1 ( Table 4). The p.(Pro290Ser) mutation is located in the sixth transmembrane (TM) domain. TM domains are important for ligand binding and receptor trafficking (Tan et al., 2004). Functional assays predicted that this genetic variation impairs not only cell surface-targeting of the receptor but also its ability to activate the Gq Protein (Monnier et al., 2009). In the present study, p.(Pro290Ser) in the PROKR2 gene was identified both in the proband and her healthy father and brother, all at homozygous states. The identification of homozygous deleterious genetic variations among healthy relatives has been previously described for diverse genetic diseases (Bonyadi et al., 2009;Biegstraaten et al., 2011;Suk et al., 2011). The question that arises is how to explain that such a deleterious genetic variation could be seen in unaffected relatives. The most plausible explanations are (i) possible digenic or oligogenic interactions due to the cooccurrence of additional genetic variants (Sykiotis et al., 2010) and (ii) the involvement of epigenetic factors (Abreu et al., 2010;Dode and Hardelin, 2010;Noel and Kaiser, 2011;McCabe et al., 2013).
Several studies employing next-generation sequencing have revealed that more and more variations known to be causative of genetic diseases can also be found among controls, thus further demonstrating the complexity of heredity (Olson, 2008;Klassen et al., 2011;Nishiguchi and Rivolta, 2012). Hence, the correlation between oligogenicity and the presence of a severe phenotype is often unclear. Sykiotis et al. (2010) reported that the same CHH phenotype seen in a propositus, carrying digenic genetic variations in FGFR1 and PROKR2 genes, was also observed in another family member who harbored only the genetic variation in the FGFR1 gene. This finding proves that carrying more or fewer genetic variations may not always correlate with the severity of disease expression, which further demonstrates the incomplete penetrance of many suspected variants. Following exome sequencing in the family of HH1, we identified 42 additional variants in 25 candidate genes in the index case as well as 39 and 33 variants in her asymptomatic brother and father, respectively. The disease-causing digenic combination profiles between the HH1 family members were compared (Figure 4). The results showed that, in addition to the PROKR2 gene variant [p.(Pro290Ser)], a second variant in c.2803C > T in the CCDC141 gene was involved in the second highest number of pathogenic digenic combinations (15%), with 18 other variants in 13 genes. The CCDC141 variant was found at a homozygous state in the patient HH1 and at a heterozygous state in the asymptomatic cases. Our analysis indicated that the zygosity state of the c.2803C > T variant in the CCDC141 gene considerably influenced the rate of pathogenic combinations. The same variant, when heterozygous in the asymptomatic cases, contributed to only 3% of the total number of pathogenic digenic combinations (Figure 4). In the same context, the implication of the IL17RD and PCSK1 variants in inducing the expression of CHH was ruled out in our propositus as they took part in the same number of digenic combinations with similar classification scores in all three family cases (Figure 4). The CCDC141 gene encodes for a coiled-coil domain-containing protein whose function is not clearly established. In mouse models, it has been shown that ccdc141 is expressed in GnRH neurons and olfactory fibers. The role of CCDC141 in the embryonic migration of GnRH neurons has been highlighted in human patients and mouse models (Hutchins et al., 2016). The prevalence of CCDC141 variants in HH patients was estimated to be 3.3% in the cohort studied by Turan et al. (2017). A nonsense variant p.(Arg724 * ) was identified in 20 probands with KS along with a loss-offunction variant in the FEZF1 gene (Hutchins et al., 2016). In 75% of the reported families, at least one additional likely pathogenic mutation in another causative gene was identified, hence underscoring the oligogenic involvement of CCDC141 gene variants in CHH (Turan et al., 2017). The c.2803C > T variant identified in the HH1 family has not been previously associated with CHH. However, it has previously been annotated as damaging in genome-wide association studies (Marsman et al., 2014;van den Berg et al., 2017;Lin et al., 2018). In light of the in silico digenic effect prediction, we hypothesize that its implication in 18 disease-causing combinations, compared to only two pathogenic variant pairs in the healthy family members, may be at the origin of disease expression in the index case (HH1), while suggesting that the full penetrance of the CCDC141 gene variant is associated to an autosomal recessive inheritance mode. It is also worth noting that the brain structure and function present differences between males and females. More than 2,000 genes show expression divergence between the two sexes at all developmental stages, especially during puberty (Shi et al., 2016). This may suggest that the variable expressivity in the HH1 family might be, in part, related to the different sex-specific gene expression profiles between our female patient and the asymptomatic male family members. For the HH12 family, the healthy mother and her affected son were both heterozygous for the p.(Lys205del) variation in the PROKR2 gene. The same situation was previously described among several patients (Dodé et al., 2006;Abreu et al., 2008;Cole et al., 2008;Sykiotis et al., 2010). The WES analysis revealed a total of 33 variants [33 variants in coding regions/splice sites (three rare); five in UTR] in 28 genes. The ORVAL prediction revealed five disease-causing digenic combinations involving DUSP6, ANOS1, DCC, PLXNA1, PROP1, and SEMA7A genes (  (Glu436Lys)] was predicted as VUS by Varsome. Sanger validation revealed the absence of this mutation in the healthy mother. The SEMA7A and DUSP6 genes were implicated in a digenic combination classified as "dual molecular diagnosis" by ORVAL. We hypothesize that the disease expression in HH12 may be explained by the digenic transmission of the SEMA7A and DUSP6 variants or a monogenic inheritance involving only the SEMA7A VUS if further functional assays would allow its reclassification into pathogenic.

CONCLUSION
The presence of homozygous PROKR2 deleterious genetic variations in asymptomatic first-degree relatives and siblings is apparently insufficient to cause KS/nIHH. This strongly suggests the incomplete penetrance trait of the disease. Further investigations are required to clarify the involvement of additional genetic and environmental factors in GnRH deficiency. Exome data analysis allowed us to explore the implication of oligogenism as a possible mechanism involved in the incomplete penetrance witnessed in the investigated families. Nonetheless, the molecular mechanisms that modulate the oligogenic interactions are far from being elucidated since the exact roles of some susceptibility genes in the physiological regulation of the GnRH nervous system are yet to be discovered.

DATA AVAILABILITY STATEMENT
Processed data related to Sanger sequencing and digenic combinations are available in the article. Raw data related to WES are available from the corresponding author upon reasonable request. Note that in Tunisia, genetic data are considered as personal private data, therefore the minimal dataset was submitted as Supplementary Material. The complete raw data may be made available upon request after obtaining IRB approval.

ETHICS STATEMENT
Patients were interviewed by both a clinical investigator and a geneticist using a structured hypogonadotropic hypogonadism questionnaire. Written informed consent was obtained from both patients and controls in accordance with the Declaration of Helsinki and after being approved by the Institutional Review Board (Registration numbers IRB00005445 and FWA00010074).

AUTHOR CONTRIBUTIONS
RM contributed to exome analysis, manuscript drafting, molecular investigation, quality control, and digenic score calculation. LCBA contributed to data analysis, manuscript drafting, recruitment of patients, molecular diagnosis, technical supervision, and quality control. SL contributed to functional in silico annotation and exome result interpretation. ZT and MS contributed to clinical diagnosis and characterization, and result evaluation. SE and CN contributed to whole exome data analysis and quality control of bioinformatic pipeline. YB contributed to molecular study, technical assistance, and 3D protein structure determination and interpretation. KM and AB contributed to data evaluation and quality control. SA contributed to study design and supervision, quality control, and manuscript revision. OM contributed to study design and supervision, exome data interpretation, quality control, and manuscript drafting and revision. All authors contributed to the article and approved the submitted version.

FUNDING
This work has been funded by the Tunisian Ministry of Public Health, the Ministry of Higher Education and Scientific Research (LR20IPT05), and the E.C. grant agreement no. 295097 for FP7 project GM-NCD-Inco.

ACKNOWLEDGMENTS
The authors would like to thank the patients for their cooperation throughout this study. The authors would also like to thank Catherine Dodé and Suhair Sunoqrot for their fruitful comments.