Broad Phenotypes of Disorders/Differences of Sex Development in MAMLD1 Patients Through Oligogenic Disease

Disorders/differences of sex development (DSD) are the result of a discordance between chromosomal, gonadal, and genital sex. DSD may be due to mutations in any of the genes involved in sex determination and development in general, as well as gonadal and/or genital development specifically. MAMLD1 is one of the recognized DSD genes. However, its role is controversial as some MAMLD1 variants are present in normal individuals, several MAMLD1 mutations have wild-type activity in functional studies, and the Mamld1-knockout male mouse presents with normal genitalia and reproduction. We previously tested nine MAMLD1 variants detected in nine 46,XY DSD patients with broad phenotypes for their functional activity, but none of the mutants, except truncated L210X, had diminished transcriptional activity on known target promoters CYP17A1 and HES3. In addition, protein expression of MAMLD1 variants was similar to wild-type, except for the truncated L210X. We hypothesized that MAMLD1 variants may not be sufficient to explain the phenotype in 46,XY DSD individuals, and that further genetic studies should be performed to search for additional hits explaining the broad phenotypes. We therefore performed whole exome sequencing (WES) in seven of these 46,XY patients with DSD and in one 46,XX patient with ovarian insufficiency, who all carried MAMLD1 variants. WES data were filtered by an algorithm including disease-tailored lists of MAMLD1-related and DSD-related genes. Fifty-five potentially deleterious variants in 41 genes were identified; 16/55 variants were reported in genes in association with hypospadias, 8/55 with cryptorchidism, 5/55 with micropenis, and 13/55 were described in relation with female sex development. Patients carried 1-16 variants in 1-16 genes together with their MAMLD1 variation. Network analysis of the identified genes revealed that 23 genes presented gene/protein interactions with MAMLD1. Thus, our study shows that the broad phenotypes of individual DSD might involve multiple genetic variations contributing towards the complex network of sexual development.


INTRODUCTION
Disorders/differences of sex development (DSD) occur when there is a discordance between chromosomal, gonadal, and genital sex (Ostrer, 2014). DSD may be due to mutations in any of the genes involved in sex determination and development in general, as well as gonadal and/or genital development specifically (Ostrer, 2014).
However, the role of MAMLD1 in sex development is controversial for several reasons: a) some MAMLD1 variants are present in the normal population (Fukami et al., 2006;Chen et al., 2010;Gaspari et al., 2011;Kalfa et al., 2011); b) the same MAMLD1 variant may be present in patients with different phenotypes (Camats et al., 2015); c) MAMLD1 variants are not present in all DSD individuals of the same family (Fukami et al., 2006); d) several MAMLD1 mutations present wild-type activity in functional studies (Camats et al., 2015); and e) the Mamld1-knockout male mouse presents with normal genitalia and reproduction (Miyado et al., 2012;Miyado et al., 2017).
MAMLD1 is expressed in human fetal and adult testis and in human ovaries (Fukami et al., 2006;O'Shaughnessy et al., 2007;Camats et al., 2015), and seems to be involved in sex development in fetal life and in adult reproductive function. Yet its exact role is not clear. It is expressed in mice gonadal cells during start of androgen biosynthesis up to male external genitalia formation and is therefore thought to be involved in the expression of Leydig-cell genes (Miyado et al., 2012), as well as supporting testosterone production in critical periods of male development (Fukami et al., 2008;Nakamura et al., 2011). In contrast, Mamld1-KO mice present normal external genitalia (but small testes and reduced seminiferous tubule size and proliferating germ cells) and reproduce similarly to wild-type mice (Miyado et al., 2012;Miyado et al., 2017). These findings challenge the role of MAMLD1 in sex development.
In a previous study, we tested functional activity of nine MAMLD1 variants detected in nine 46,XY DSD patients with broad phenotypes (Camats et al., 2015). None of the MAMLD1 mutants, except truncated L210X, had diminished transcriptional activity on known target promoters CYP17A1 and HES3. In addition, protein expression of MAMLD1 variants was similar to wild-type, except for the truncated L210X. We therefore hypothesized that MAMLD1 variants may not be sufficient to explain the phenotype in 46,XY DSD carriers, and that further genetic studies should be performed to search for additional hits explaining the broad variability.
In the past decade, High throughput sequencing (HTS) has changed the genetic approach in research and diagnostics. Whole-exome sequencing (WES) has led to the discovery of many new genes and has given insight into complex traits. Oligogenic inheritance is currently discovered for several disorders by HTS. In the field of sex development, digenic inheritance has recently been suggested in a 46,XY DSD patient with gonadal dysgenesis (NR5A1 and MAP3K1 variants) (Mazen et al., 2016); in a family with 46,XY DSD males (NR5A1 variants) and 46,XX POF females (NR5A1 and TBX2) (Werner et al., 2017); as well as in a DSD patient with ambiguous genitalia, micropenis, and inguinal testes (SEMA3A and AKR1C4) (Fan et al., 2017). Similarly, we found oligogenic origin of disease in heterozygous NR5A1 46,XY DSD patients by performing WES (Camats et al., 2018). In addition, in patients with hypospadias, an oligogenic origin was suggested by two other NGS studies (Kon et al., 2015;Eggers et al., 2016).
Therefore, in this study, we performed WES in seven 46,XY patients with DSD (Camats et al., 2015) and one 46,XX patient with ovarian insufficiency, who all carried MAMLD1 variants. WES data were filtered by common tools and a disease-tailored algorithm including MAMLD1-related and DSD-related known and candidate genes. Additional hits in likely diseasecausing genes were detected in all eight MAMLD1 carriers. Our results suggest that oligogenic origin of disease may contribute towards the broad phenotypes of human MAMLD1.

Patients
The study was approved by the Ethics Committee of Hospital Universitari Vall d'Hebron (Barcelona, Spain) (CEIC: PR(IR)23/2016). Written informed consent was obtained from the patients for the publication of their cases. Eight DSD patients (seven 46,XY and one 46,XX) each carrying one MAMLD1 variant were analyzed using WES. Clinical and genetic characteristics of 46,XY patients were previously reported in detail in Camats et al., (2015) and are summarized in Table 1 together with the 46,XX patient.

DNA Extraction, WES and Bioinformatic Analysis
DNA was extracted from blood leukocytes using QiaCube (Qiagen, Hilden, Germany) or manually using a DNA isolation kit (Qiagen). WES was performed by CNAG (Centre Nacional d' Anàlisi Genòmica, Barcelona, Spain). Libraries were prepared with a SureSelect Human All Exon V5 capture kit (Agilent, Santa Clara, CA, USA) and sequenced with a HiSeq™ 2000 sequencing system (v3, 2x100, Illumina, San Diego, CA, USA). Putative candidate variants were confirmed by Sanger sequencing. The genomic datasets were annotated (alignment with human genome hg19/grch37) and filtered with the functional annotation of genetic variants from HTS data (ANNOVAR; http://annovar. openbioinformatics.org/) (Wang et al., 2010), visualized and explored in Integrative Genomics Viewer (IGV, Broad Institute, Cambridge, MA, USA; https://www.broadinstitute.org/igv/ (Robinson et al., 2011). Frequencies of variants of relevant candidate genes were obtained from the Genome Aggregation Database (gnomAD; http://gnomad.broadinstitute.org/about) (Lek et al., 2016) and the Collaborative Spanish Variant Server (CSVS; CIBERER BIER, Valencia, Spain; http://csvs.babelomics. org/; August 2018) (Dopazo et al., 2016). gnomAD includes gene variants from exome and genome sequencing data: 123,136 exomes and 15,496 genomes from unrelated individuals (from population and disease-specific studies). CSVS database includes (among others) exomes from a population of 267 healthy unrelated subjects. WES data were filtered by a disease-tailored list of MAMLD1related and DSD-related known and candidate genes (n = 606) similar to the algorithm previously set up for Camats et al., (2018). We generated a project-specific filter for DSD-related and MAMLD1-related genes by searching in published literature and databases. DSD-related genes are part of our DSD-gene database and tools (Camats et al., 2018), which have been currently updated. The DSD-related gene list included genes with reported (potentially) deleterious variants in patients with 46,XY and 46,XX DSD, genes with reported (potentially) disease-causing variants in syndromic patients with involvement of sex development, those "related" to DSD conditions in KO/ mutant animal models (mice and rats), and also overexpressed, upregulated or downregulated genes in rodent embryonic gonadal cells (Camats et al., 2018). For the search for functional human partners of MAMLD1 and for possible interactions within interesting genes, the Search Tool for the Retrieval of Interacting Genes/Proteins (STRING, http://string-db.org/) (Jensen et al., 2009) and the Biological General Repository for Interaction Datasets (BioGRID, thebiogrid.org) (Stark, 2006) were used.

Variant Analysis Per Patient
After annotation, variant analysis was performed by the following steps. A) Each patient's WES data were first filtered by our MAMLD1-and DSD-related known and candidate genes. B) We kept variants with MAF (minor allele frequency) ≤0.015 or not detected in gnomAD, and variants with the following predicted type, consequences and locations: splicing (intronic or exonic), exonic, intergenic, regulatory. C) We confirmed the correct annotation and location of variants by checking their alignment data in IGV (alignment with human genome hg19/ grch37) (data not shown). D) We excluded variants that were considered non-relevant for our study: E.g. 1) variants found in more than two patients, 2) variants in repeat regions, 3) variants in genes or gene regions with high variability, 4) variants with low coverage and/or low quality, 5) variants with non-similar allelic depths. E) We revised variants with the annotated pathogenic predictors: functional exonic, evolutionary-conservation and splicing predictors (ANNOVAR and Alamut Visual software), as previously described. F) We run InterVar and VarSome to classify the variants, searched for reported (potentially) human diseasecausing variants with the HGMD, and revised evidences of relationship with DSD, sex development and clinical phenotype of each patient with literature and database search. G) We used STRING to find out interactions among genes carriers of interesting variants (DSD-related and/or MAMLD1-related) (Figures 1 and 2). H) We checked MAF in a healthy cohort of Spanish population (CSVS: 267 unrelated healthy controls). I) We rejected variants with MAF ≥ 0.01 (gnomAD, CVSV, August 2018), thus less plausible to be a DSD-causing variant. Importantly, synonymous variants were not rejected because it has been shown that they may affect splicing.

RESULTS
WES performed in eight unrelated subjects (seven 46,XY and one 46,XX) harboring hemizygous/heterozygous MAMDL1 variants revealed several candidate gene variants that potentially contribute to each patient's phenotype. A detailed summary of patients' characteristics and number of variants and genes is shown in Table 1. A list of identified candidate variants and corresponding information from literature is given in Table 2.
We identified a total of 55 potentially deleterious/candidate heterozygous/hemizygous variants in 41 genes in the eight hemizygous/heterozygous MAMLD1 patients (Tables 1 and 2). In the seven 46,XY patients 1-16 variants were found in a total of 1-16 genes, while the 46,XX MAMLD1 patient revealed 14 additional variants in 13 genes. (Tables 1 and 2).
Patient 1 harbored nine variants in seven genes: CYP1A1, EVC, GRID1, NOTCH1, RET, RIPK4 and ZBTB16, all of them associated to gonadal/genital anomalies ( Table 2). Patient 2 carried one variant in RECQL4, associated with syndromic hypospadias ( Table 2). Patient 3 presented two variants in two genes: GLI2 (associated to gonadal/genital anomalies) and RECQL4 (associated to syndromic hypospadias) ( Table 2). Patient 4 had four variants in four genes: CDH23, COL9A3, MAML1 and NOTCH1; all, except MAML1, have been proposed to be associated with gonadal development ( Table 2). In patient 5, six variants in six genes were found: BNC2, FGF10, HSD3B2, IRX5, MAML2 and NOTCH2; all, except MAML2, have also been associated with hypospadias or gonadal development ( Table 2). Patient 6 carried 16 variants in 16 genes: ATF3, BNC2, CYP1A1, EYA1, FLNA, FRAS1, GLI3, HOXA13, IRX5, IRX6, MAML1, NRP1, MAML3, PROP1, PTPN11 and WDR11 ( Table 2). Thirteen of these genes are associated with risk of hypospadias and/or syndromes that include abnormal gonadal/genital development, whereas MAML1 is unrelated, MAML3 has been proposed to be associated with female gonadal development and PROP1 has only been associated with anterior pituitary insufficiency/ hypogonadotropic hypogonadism. In addition, six of these genes have previously been described in patients with aortic diseases and cardiopathies (   Nodes represent proteins. Filled nodes show proteins with known or predicted 3D structure. Empty nodes depict proteins with unknown 3D structure. Candidate genes are underlined. Known interactions correspond to curated databases (turquoise lines) and experimentally determined interactions (pink lines). Predicted interactions correspond to gene neighborhood (green lines), gene fusions (red lines) and gene co-occurrence (blue lines). Other interactions correspond to text mining (yellow lines), co-expression (black lines) and protein homology (violet lines). Genes with no interactions are on the right side of each network.

7
Month 2019 | Volume 10 | Article 746 Frontiers in Genetics | www.frontiersin.org       Among them, only MAML2 has not been related to gonadal or genitourinary system development ( Table 2).
We performed interactome analysis for the identified DSD genes using bioinformatic tools for the analysis of possible geneprotein interactions. The network comprising all genes identified is shown in Figure 1. Overall, a connection was found for 27 of the 41 genes. MAMLD1 connects directly to MAML1/2/3. Via NOTCH1/2 8 genes are in connection with MAMLD1, namely WNT9A/9B, GLI2/3, FGF10, RET, PROP1 and NRP1. Some of these genes are also central nodes for further connections; e.g. GLI3 for EVC, FGF10, GLI2, RIPK4 and EYA1; and RET for PIK3R3 with PTPN11, which also is connected with RIPK4. RIPK4 itself is a central node for ZBTB16, CUL4B, GLI3 and PTPN11. NRP1 is connected to FLNA and EYA1 connects with FRAS1 and FREM2. In addition, 2 isolated gene couples have been revealed by our analysis: CYP1A1-HSD3B2 and MYO7A-CDH23. These observations give an idea of the complex interactions among genes related to sex development.

DISCUSSION
Sex development is a very complex biological event which requires the concerted collaboration of a large network of genes in a spatial and temporal correct fashion. In the past, much has been learned about human sex development from monogenic DSD, but the broad spectrum of phenotypes in numerous DSD individuals remains a conundrum. Oligogenic disease has been proposed. In fact, multiple genetic hits, which might not be deleterious by themselves, have been found in several individuals with DSD (Kon et al., 2015;Eggers et al., 2016;Mazen et al., 2016;Werner et al., 2017;Camats et al., 2018). In a previous study of 46,XY DSD patients carrying MAMLD1 variants, we showed that none of the variants were functionally pathogenic except for a stop variant (Camats et al., 2015). In the present study, we searched for additional genetic hits in DSD patients harboring MAMLD1 mutations and manifesting with unexplained broad phenotypes. Using HTS and a custom-made  (Jameson et al., 2012, Clement et al., 2007, Beverdam and Koopman, 2006 8 Hauser syndrome (Waschk et al., 2016); bicornuate uterus (Waschk et al., 2016), organogenesis urogenital system (Carroll et al., 2005;Grinspon and Rey, 2014) All patients were heterozygote for these variants and were checked and confirmed through IGV software ( The seven MAMLD1 patients with 46,XY DSD presented phenotypes from female external genitalia (patient 4) to variable degrees of hypospadias, cryptorchidism and small penis (Table 1). Interestingly, patient 4 with female external genitalia had normal T secretion. Similarly, patient 5 carrying a heterozygous HSD3B2 variant, had normal levels of 17OH-pregnenolone, DHEA and DHEA-S (data not shown). Patient 6, who presented with a right aortic arch, was found to carry variants in five genes (BNC2, FLNA, MAML1, NRP1 and PTPN11) that have been previously described in patients with heart and/or vascular anomalies (Tartaglia et al., 2002;Bhoj et al., 2011;Lee et al., 2014;Shaheen et al., 2015;Preuss et al., 2016;Chen et al., 2018). The 46,XX patient (patient 8), with primary amenorrhea, hypergonadotropic hypogonadism, normal female external genitalia and small uterus harbored gene variants involved in gonadal development and DSD (CUL4B, DAPK1, EMX2, FREM2, IGFBP2, MAML3, MYO7A, NOTCH1, PIK3R3, TGFBI, WNT9A and WNT9B; Table 2). Five of these genes (DAPK1, IGFBP2, MAML3, PIK3R3 and WNT9A) have so far only been related to female gonadal development (Table 2).
Overall, the genes detected in our eight studied patients with MAMLD1 variants have been previously reported in humans with hypospadias, cryptorchidism, micropenis, and other urogenital abnormalities; or they have been found involved in sexual and gonadal development. Also, some of them have been associated with specific syndromes in patients with genitourinary anomalies: CAKUT syndrome, Ellis-van Creveld syndrome, Fraser syndrome 1, Fraser syndrome 2, hand-foot-genital/ Guttmacher syndrome, Noonan syndrome, Mayer-Rokitansky-Küster-Hauser syndrome, Popliteal pterygium syndrome and Rothmund-Thomson syndrome (see Table 2). However, none of the present patients presented a complete phenotype for any of these syndromes, maybe because none of the variants completely impairs gene expression and protein function, as inferred by the in silico analyses. Detailed information on these genes from current literature is given in Supplementary Materials (S1).
A search for an underlying network comprising variants in the identified genes related to MAMLD1 revealed a considerable number of genes which showed gene-gene, gene-protein or protein-protein interactions (Figures 1 and 2) suggesting that genetic variations in these genes may affect sex development. In addition, MAML3 was found in a network related to female gonadal development (Jameson et al., 2012). Accordingly, one variant in MAML3 was present in our 46,XX patient. The analysis of gene/protein network interactions per patient gives an idea of the complexity of the interactions among genes related to sex development. The more variants detected in DSD-related genes, the better to build an interaction network searching for clues on genetic relationship(s) for sex development. In our DSD individuals carrying MAMLD1 variants, three genes seemed prominent in the network analysis, NOTCH1, NOTCH2 and GLI3. NOTCH signaling is a highly conserved signaling pathway and comprises 4 transmembrane receptors. It is essential for the regulation of embryonic development of multiple organ systems including gonadal development (Windley and Wilhelm, 2016). NOTCH signaling is implicated in Leydig cell differentiation in an inhibitory regulatory fashion (Windley and Wilhelm, 2016). Autosomal dominant mutations in NOTCH1 cause the Adams-Oliver syndrome (OMIM 616028), while autosomal dominant mutations in NOTCH2 are reported in the Alagille syndrome 2 (OMIM 610205) and in the Hajdu-Cheney syndrome (OMIM 102500). By contrast, GLI3 is a zinc-finger transcription factor belonging to the desert hedgehog (DHH) signal transduction pathway. DHH signaling is essential for driving Leydig cell differentiation (Windley and Wilhelm, 2016). Thus, NOTCH and DHH signaling work together to regulate Leydig cell development (Windley and Wilhelm, 2016). Autosomal dominant mutations in GLI3 are described in the Pallister-Hall syndrome (OMIM 146510) or in the Greig cephalopolysyndactyly syndrome (OMIM 175700).
Taken together, our results expand the landscape of genes possibly involved in DSD by revealing both new and old players. Genetic platforms for DSD diagnostics currently consider about 270 genes that have been identified with monogenetic forms of DSD in (mostly) several independent individuals (Cools et al., 2018). Our eight MAMLD1 individuals share variants in 19 genes comprised in such DSD panels, including ATF3, BNC2,CUL4B,EVC,FLNA,FRAS1,FREM2,GLI3,HOXA13,HSD3B2,IRX5,NOTCH2,PROP1,PTPN11,RECQL4,RET,RIPK4,WDR11 and ZBTB16. By contrast, through our work 22 new genes are now added for considering with differences in sex development: CDH23, COL9A3,CYP1A1,DAPK1,EMX2,EYA1,FGF10,GLI2,GRID1,IGFBP2,IRX6,MAML1,MAML2,MAML3,MYO7A,NOTCH1,NRP1,PIK3R3,PPARGC1B,TGFBI,WNT9A and WNT9B. Ideally, genetic variants are tested functionally for proof of their disease-causing effect in model systems. However, when finding multiple variants, which may all contribute only partially, such testing is no longer feasible. Therefore, the likelihood of disease-causing effect of identified variants was assessed in our study by established bioinformatic tools for genetics and by assessing the genotype-phenotype correlation in each patient with current knowledge from literature and databases in the field. In future studies with bigger sample size, next-generation statistical genetic analyses may be employed to identify associations between a group of variants and the complex trait of sex development (Weissenkampen et al., 2019).
In summary, HTS analysis indicates that the broad DSD phenotypes of MAMLD1 patients may be due to additional variants in other DSD-related genes. We found up to 55 additional genetic hits that may contribute to the DSD phenotype making an oligogenic causation plausible. Bioinformatic network analysis can help in interpreting complex genetic data and put identified single candidate genes into a greater perspective to understand their possible role in DSD biology.

ETHICS STATEMENT
The study was approved by the Ethics Committee of Hospital Universitari Vall d'Hebron (Barcelona, Spain) (CEIC: PR(IR)23/2016). Written informed consent was obtained from the patients for the publication of their cases.

AUTHOR CONTRIBUTIONS
CF: Conceptualization, funding acquisition, investigation, methodology, interpretation, project administration, resources, supervision, writing -original draft preparation, writing -review and editing. LA: Interpretation, supervision, resources, visualization, writing -review and editing. MF-C: Investigation, writing -review and editing. K-SS: Validation, interpretation, writing -review and editing. IM: Resources, interpretation, writing -review and editing. LC: Resources, writing -review and editing. IE: Resources, writing -review and editing. NC: Conceptualization, data-curation, formal analysis, investigation, methodology, interpretation, project administration, supervision, visualization, writing -original draft preparation, writing -review and editing.