This study was approved by the Ethics Committee of The First Affiliated Hospital, Zhejiang University School of Medicine, China (approval number 2018–727). Written and informed consent was obtained from the parents of the proband.

The proband and her elder sister were assessed for medical history and clinical manifestations. Physical examination, especially on the external genitalia, was also performed. Thereafter, laboratory examination including hormonal analysis (follicle-stimulating hormone, luteinizing hormone, estradiol, progesterone, prolactin and testosterone), karyotype analysis, imagological examination, laparoscopy and gonadal biopsy with histopathological examination were performed.

Karyotype analysis was carried out using standard techniques. Briefly, genomic DNA was obtained from both the proband and her parents via peripheral blood leukocytes using standard procedures. We then performed whole-exon sequencing (WES) as previously described (Chen et al., 2019). The Sanger sequencing variant detected in the MAP3K1 gene was described according to the NCBI entry NG_031884.1 (NM_005921.2). The variants were named according to the Human Genome Variation Society (HGVS) sequence variant nomenclature (den Dunnen et al., 2016).

We downloaded various MAP3K1 protein sequences from Uniprot (https://www.uniprot.org/), while ClustalX program (http://ftp-igbmc.u-strasbg.fr/pub/ClustalX/) was used for sequence alignment. Sequence alignment results were displayed online using ConSurf (https://consurf.tau.ac.il/) (Ashkenazy et al., 2010; Ashkenazy et al., 2016).

The variant’s pathogenicity was analysed using the webserver; PREDICT-SNP (https://loschmidt.chemi.muni.cz/predictsnp1/) (Bendl et al., 2014). This tool employs a consensus of five different predictors, which include PredictSNP, PolyPhen-1 (Ramensky et al., 2002), PolyPhen-2 (Adzhubei et al., 2010), SIFT (Ng and Henikoff, 2003), and SNAP (Ramensky et al., 2002) along with provision of a confidence score for the predictions.

Due to lack of experimental structure or sufficient homology models, the 3D structure of MAP3K1 protein was generated via ab initio modeling. The amino acid sequence of wild-type human MAP3K1 protein was retrieved from UniProt database (http://www.uniprot.org). The 3D structure of the wild-type human MAP3K1 protein was predicted through Iterative Threading ASSEmbly Refinement (I-TASSER) (https://zhanglab.ccmb.med.umich.edu/I-TASSER/) (Roy et al., 2010; Yang et al., 2015). Thereafter, the online server Dynamut (http://biosig.unimelb.edu.au/dynamut/) was used to assess local conformation of wild-type and variant proteins, and assess changes, if any, in protein flexibility (Rodrigues et al., 2018). A structural representation was generated using PyMOL 2.4 (https://pymol.org/2/), a molecular visualization system in the open-source foundation.

Wild-type full-length human MAP3K1 complementary DNA (cDNA) (NM_005921.2) was chemically synthesized and cloned into a pCDH vector. The c.556A > G variation was introduced into the MAP3K1 sequence using an In-Fusion HD Cloning kit (TaKaRa, United States). Forward and reverse primers used for PCR amplification were 5′ CGT​CCA​GAG​GAA​CGA​ATG​ATC​GGG​GAG​AAA​CTG​AAG​GCA​A3′ and 5′ GAT​CAT​TCG​TTC​CTC​TGG​ACG​ATC​ATC​CA3’, respectively. We then constructed plasmids transiently expressing MAP3K1-3×Flag, MAP3K1-3×HA, RHOA-3×Flag, or UB-3×Ha were for co-immunoprecipitation assays. The MAP3K1 gene was cloned into the plasmid p3×Flag-CMV-14 or pcDNA3.1(+)-3×HA plasmid. RHOA (NM_001664) gene was cloned into the p3×Flag-CMV-14 plasmid, while UB (NM_021009.7) gene was cloned into the pcDNA3.1(+)-3×HA plasmid. The ligated product was transformed into Escherichia coli (DH5α cells). The recombinant plasmid clones were sequenced to verify the integrity of the insertion sequences.

HEK-293T cells were co-transfected with TCF/LEF-TOP flash reporter plasmid and wild-type or variant MAP3K1 gene expression vectors. The cell transfection experiments were performed in a 24-well plate. The plasmids were transfected into HEK-293T cells using Lipofectamine 3,000 (Thermo Fisher, United States), according to the manufacturer’s instructions. After each variant was co-transfected, the DNA ratio of the MAP3K1 and TCF/LEF-TOP flash reporter plasmids (Genechem, China) were 1:1 (250 ng each). After transfection, we employed the luciferase Reporter Assay System (Promega, United States) at 24 h, 36 h, 48 h, 60 h, and 72 h. The transfection efficiency was quantitatively assessed using quantitative PCR (qPCR), as shown in Supplementary Figure S4. The experiments were independently repeated at least five times, with at least three replicates. Data were expressed as mean ± SEM and analyzed using a t-test.

The structural complementation reporter system, NanoLuc^® Binary Technology (NanoBiT) (Promega, Cat#N2014), was used to analyze protein interactions. The NanoBiT consists of Large BiT (LgBiT; 18 kDa) subunit and a small complimentary peptide (SmBiT; 3.6 kDa) (Dixon et al., 2016). When the two target proteins interact, the LgBiT and SmBiT subunits converge to form an active enzyme, which produces a bright luminescent signal in the presence of a substrate. For NanoBiT protein interaction experiments, the MAP3K1 gene CDs and the N-terminus of the LgBiT or the C-terminus of the SmBiT were fused to form MAP3K1-LgBiT-N or MAP3K1-SmBiT-C plasmids, respectively. The RHOA gene and the N-terminus of the LgBiT were fused to form RHOA-LgBiT-N expression plasmid. On the other hand, the UB gene and the C-terminus of the LgBiT were fused to form the UB-LgBiT-C expression plasmid. HEK293T cells were then co-transfected with MAP3K1-SmBiT-C and RHOA-LgBiT-N to assess the ability of the MAP3K1-WT or MAP3K1-Mut to interact with RhoA. Similarly, the HEK293T cells were co-transfected with SmBiT-MAP3K1-N and LgBiT-UB-C to evaluate the ability of the MAP3K1-WT or MAP3K1-Mut to interact with Ub. The HEK-293T cells in a 24-well plate were transfected with a total of 500 ng DNA using a Lipofectamine™ 3,000 transfection reagent. Twenty-4 hours after transfection, the culture medium was replaced with fresh OPTI-MEM, and then a diluted substrate was directly added to each well. Fluorescence intensity was measured using a Multimode Plate Reader VICTORNivo (ND-1000, Thermo Fisher). The transfection efficiency was quantitatively analyzed through qPCR, as shown in Supplementary Figure S5. The experiments were independently repeated at least three times, with at least three replicates. Data were expressed as mean ± SEM and analyzed using a t-test.

Co-IP was performed as previously described (Wang et al., 2019). Briefly, the MAP3K1-3×HA was co-transfected with RHOA-3×Flag, while MAP3K1-3×Flag was co-transfected with UB-3×HA plasmid using Lipofectamine 2000. At 48h post-transfection, the cells were collected and washed with PBS, and then harvested in ice-cold modified RIPA lysis buffer. Cell lysates were immunoprecipitated, and then incubated overnight with anti-Flag or anti-HA magnetic beads (GenScript, China) at 4°C. The proteins were denatured at 95°C for 10 min and then resolved in SDS-PAGE. Thereafter, the blots underwent western blot analysis with anti-HA or anti-FLAG antibodies. Protein band intensities were determined and used to calculate the relative ratio of Co-IP to IP signal. Data were shown as mean ± SEM (n = 3) and analyzed using Mann-Whitney U test.

The HEK293T cells and a human pluripotent teratocarcinoma testicular origin malignant embryonal carcinoma cell line (NT2/D1, ATCC^® CRL-1973™) were cultured in HG-DMEM (Invitrogen, United States), containing 10% fetal bovine serum (FBS). On the other hand, the ovarian origin granulosa cell tumor line (KGN, RIKEN Cell Bank^® RCB1154) was cultured in DMEM/F12 medium (Invitrogen, United States), containing 10% FBS. The cells were grown at 37°C with 5% CO₂.

Total RNA was obtained using Total RNA Extraction Reagent (Promega). Thereafter, mRNA expression levels were analyzed by real-time PCR using the ABI PRISM 7900HT Sequence Detection System (Applied Biosystems) (Liu et al., 2014). Expression values were normalized using β-actin (ACTB gene). Supplementary Table S1 shows the basic information on primers, gene name, gene ID and NCBI reference sequence.

The cells were collected and lysed in RIPA lysis buffer. Equal amount of protein (25 ug), as quantitated using the BCA Protein Assay, from each group was resolved in 10% SDS-PAGE and then electro-transferred onto a polyvinylidene difluoride membranes for western blott analysis (Li et al., 2018). The sources and catalog numbers of the antibodies for Western blotting are as listed in Supplementary Table S2. Each experiment was independently performed at least three times.

Western blot data were processed using the Fiji/ImageJ software. Statistical analyses were performed in Graph Pad Prism 7 program. Statistical significance was tested using a t-test or a Mann-Whitney U test. The results were drawn from at least three independent experiments. Data were shown as mean ± SD (n ≥ 3). * indicates p < 0.05, ** indicate p < 0.01 and *** indicate p < 0.001.

A 13.5-year-old girl was admitted to the First Affiliated Hospital of Zhejiang University School of Medicine due to lack of breast development and menstruation. Both parents were Han Chinese, healthy and nonconsanguineous. Her height and weight were 157.8 cm (+0.03SDS) and 39 kg (-1.06SDS), respectively. Based on the Tanner scale, her breasts and pubic hair were Tanner I. External genitalia presented as a normal female with no clitoromegaly. Ultrasonography analysis revealed a rudimentary uterus (1.5 × 0.5 × 1.0 cm) and no identifiable gonadal tissue (Figures 1A,B). Analysis of the blood hormonal features of the proband demonstrated a sex steroid deficiency and high gonadotropin levels (Table 1). Further diagnostic tests demonstrated a 46,XY karyotype. Following genetic testing, the proband underwent a laparoscopic bilateral gonadectomy. In addition, stunted uterine-like tissue was observed during the operation, but no obvious abnormalities were found in the bilateral fallopian tubes. Streak gonads could be seen below the bilateral fallopian tubes. Pathological evaluation revealed that the bilateral streak gonads were ovarian interstitial-like cells and germinal epithelium, without any malignant transformation (Figure 1C).

Her elder sister, a 19-year-old female, had undergone a laparoscopic bilateral ovarian tumor removal 5 years ago. Her external genitalia presented as a normal female with no clitoromegaly. The vagina ended in a blind pouch, with normal labia. However, the vagina and urethra had separate openings. Her karyotype was also 46,XY. Ovarian ultrasonography indicated the presence of a 12.2 × 9.5 × 7.0 cm mass in the left ovary and a 2.3 × 2.64 × 1.8 cm mass in the right ovary (Figures 1D,E). Surgical records showed that she had bilateral ovarian tumors of about 20.0 × 14.0 × 7.0 cm on the left and 5.0 × 2.5 × 2.0 cm on the right. In addition, HE staining, placental alkaline phosphatase staining, and immunohistochemical staining (CD117) classified her bilateral ovarian tumor as dysgerminoma (Figure 1F, Supplementary Figure S1). Serum gonadotropin levels of the sister were elevated, and gonadal steroid levels were low, demonstrating gonadal failure. Based on the above clinical characteristics, both patients were diagnosed with 46,XY DSD.

Both WES and Sanger sequencing were performed on samples from the members of a 46,XY DSD core pedigree, including two affected siblings and their parents (Figure 2A). Quality control of FASTQ files was performed, and the proband’s WES satisfied the variant screening criteria (Supplementary Figure S2). There were a total of 21,230 variants in the exome region, out of which 53.15% were nonsynonymous variants (Figures 2B,C). The 11,282 nonsynonymous variants were filtered into 25 rare variants. Finally, a mutation site was identified in the MAP3K1 gene (c.556A > G) (Figure 2B). The MAP3K1 gene underwent amplification by PCR. The Sanger sequencing demonstrated that the proband’s mother and elder sister also carried the c.556A > G (r.556a > g, p. R186G) heterozygous variant (Figures 2A,D). Besides, the WES and Sanger sequencing confirmed a missense variant of the MAP3K1 gene, which is not contained in the public population database gnomAD (http://gnomad), Chinese Millionome database (https://db.cngb.org) or Human Gene Mutation database (http://www.hgmd.org/).

To further analyze the genetic testing data, we assessed the interspecific conservation and pathogenicity of the c.556A > G variant. We compared MAP3K1 protein sequences in 72 vertebrates such as Homo sapiens, Danio rerio, Mus musculus, Pan paniscus, Loxodonta africana, Gallus gallus, Macaca mulatta, Bos taurus, Gopherus agassizii, Xenopus laevis, Ictidomys tridecemlineatus, and Anolis carolinensi, as shown in Figure 2E. ConSurf Server was used for alignment and tint the residue according to the degree of conservation. The residues at position Arg186 of MAP3K1 were highly conserved across all the represented vertebrates (100%) (Figure 2E). The MAP3K1 missense variant was predicted to be deleterious by all five predictors (PredictSNP, polyphen1, polyphen2, SIFT, and SNAP), with a high confidence score of 0.65, 0.59, 0.60, 0.79, and 0.81, respectively (Supplementary Figure S3). These findings suggested that the c.556A > G (p.R186G) variant in MAP3K1 gene might be pathogenic.

To further understand the structural changes in the MAP3K1 protein, we used the I-TASSER online software to perform ab initio modeling of the MAP3K1 structure (Figures 3A–C). C-score, a confidence rating factor with values in [−5,2], was used to evaluate the quality of the I-TASSER predictions. Our analysis showed that the MAP3K1 model obtained by the I-TASSER server had a high C-score (0.56–0.80) and TM-score, with a structural similarity measurement with values in [0,1], correlation coefficient (0.91), thus demonstrating that the predictions were credible. The p. R186G variant was located at the middle of the GEF domain (amino acids 164–231) and in the predicted alpha-helix of the MAP3K1 protein.

The I-TASSER prediction results were then used as a basis for further modelling the protein structure to predict possible changes in protein flexibility after mutation. The DynaMut web server was used to analyze and visualize protein dynamics, and assess the impact of mutations on the protein dynamics. We showed a molecular vector field representation of the first non-trivial mode of the molecule motion based on Normal Mode Analysis, as shown in Figure 3D. Figure 3E shows deformation energy performed over the first ten non-trivial modes of the molecule. This data provides a measure for local flexibility in the wild-type MAP3K1 protein. The data showed that the flexibility of the GEF domain is highest of all domains within the MAP3K1 protein. In addition, the p. R186G variant further increased the flexibility of the GEF domain (ΔΔSVibENCoM: 0.828 kcal/mol/K) (Figure 3F).

Finally, we predicted the local interactions between the wild-type and variant amino acid residues via the DynaMut online server. The data demonstrated that the variant led to a change from arginine (Arg) to glycine (Gly), plus significant charge differences. In particular, the loss of charge of the arginine residue could potentially result in other interaction losses with other molecules or residues. In addition, there was introduction of a more hydrophobic residue into the variant, which could trigger loss of hydrogen bonds and/or disruption of protein folding (Figures 3F,G). Thus, the unique structural disturbance of the p. R186G variant altered the local hydrogen bond pattern of the MAP3K1 protein, which might be responsible for the increased flexibility of the variant.

To examine if the MAP3K1^R186G variant altered its interaction with other cofactors, we used a luciferase-based protein fragment complementation assay (NanoBiT), to measure its ability to bind RhoA or Ub. Co-expression of the SmBit-MAP3K1-C and LgBit-RhoA-N yielded a clear luminescent signal. This finding suggested a strong interaction for the receptors in this combination (Figure 4A). As shown in Figures 4B,B', the binding of the variant to RhoA was significantly increased to 3.19 ± 0.18 fold compared to the wild-type protein. The interaction between MAP3K1 and Ub was simultaneously analyzed using the same method (Figures 4C,D). The binding ability of this variant to Ub was reduced to 43–49% of that of the wild-type (Figures 4D,D'). These data indicated that MAP3K1^R186G had increased binding affinity towards RhoA but reduced affinity for Ub.

RhoA and Ub have a common binding motif. RhoA interacts with the GEF and PHD domains (residues 149–636) of MAP3K1 (Gallagher et al., 2004). On the other hand, Ub interacts with the PHD domains (residues 438–440) of MAP3K1 (Charlaftis et al., 2014). The variant site (p.R186G) is located in the RhoA and MAP3K1 binding regions. A previous study showed that a combination of MAP3K1 and RhoA may inhibit MAP3K1-specific ubiquitination of other substrates, including auto-ubiquitination (Gallagher et al., 2004). Our Co-IP assay confirmed that the binding capacity of MAP3K1^R186G to RhoA was 2.57 times that of the wild-type protein (Figures 5A,A'), and the ubiquitination degree of the wild-type MAP3K1 protein was 3.6-fold that of the MAP3K1^R186G variant (Figures 5B,B').

RhoA is a potent activator of MAP3K1 through activation of its kinase (Gallagher et al., 2004). Glycogen synthase kinase 3β (GSK3β) is a classic negative regulator of Wnt/β-catenin signaling pathway and phosphorylation of GSK3β results in its inactivation. p38 regulates canonical Wnt/β-catenin signaling by inactivation of the GSK3β (Bikkavilli et al., 2008). To determine the effect of the MAP3K1^R186G variant on kinase and ovarian-specific pathways, we evaluated phosphorylation levels of p38 and GSK3β as well as activation of the Wnt4/β-catenin pathway in the testicular teratoma cells (NT2/D1) and ovary-derived granular cells (KGN). The NT2/D1 and KGN cells were transiently transfected with expression plasmids, containing either the wild-type or mutant MAP3K1 gene. The MAP3K1 protein levels were significantly elevated in the mutant group compared to the wild-type group. The protein expression of the MAP3K1^R186G variant was 6.4 and 2.0 times higher than that of the wild type in the KGN and NT2/D1 cells, respectively (Figures 6A–D). In the KGN ovarian cells, the phosphorylated p38 (P-p38) and GSK3β (P-GSK3β) protein levels were 2.5 and 2.6-fold higher, respectively, in the mutant group compared to those in the wild-type group. Moreover, our results demonstrated that the variant could activate the Wnt4/β-catenin signaling. The mutant group showed 1.7 and 3.7-fold increase in the Wnt4 and β-catenin protein levels, respectively, compared to the wild-type group in the KGN cells (Figures 6A,B). Similarly, in the NT2/D1 testicular cells, phosphorylated p38 and GSK3β proteins were significantly elevated in the mutant group, which were 2.3 and 2.2 times higher than those in the wild-type group, respectively. Wnt4 and β-catenin levels had also increased to about 1.7 and 2.5-fold in the mutant group compared with the wild-type group (Figures 6C,D). Furthermore, we analyzed the transcriptional activity of the variant on β-catenin using a TCF/LEF luciferase assay. The data showed that at 48 h after transfection, the relative activity of TCF/LEF in the mutant group reached its peak, which was 9.78 ± 0.57 times that of the wild-type group (Figures 6E,F). These results suggested that the MAP3K1^R186G variant enhances the phosphorylation of p38 and GSK3β, leading to the accumulation of β-catenin and activation of the Wnt4/β-catenin signaling pathway.

To further assess the effect of the MAP3K1^R186G variant on sex-specific signaling pathways, we analyzed the effects of the variant on the expression of critical genes in testis and ovarian differentiation. Forkhead box L2 (FOXL2), a pro-ovarian transcription factor, is essential for ovarian differentiation and functions (Elzaiat et al., 2017; Nicol et al., 2019). Fibroblast growth factor 9 (FGF-9) and its receptor, fibroblast growth factor receptors 2 (FGFR-2), as well as doublesex and mad-3 related transcription factor 1 (DMRT1), are critical regulators of testicular determination and mediate inhibition of ovarian pathways (Bagheri-Fam et al., 2017; Huang et al., 2017; Harpelunde Poulsen et al., 2019). We overexpressed the wild-type or mutant MAP3K1 gene in the KGN ovarian granulosa cell line. Our analysis showed that there was no significant difference in MAP3K1 mRNA expression between the wild-type and mutant groups (Figures 7A,C). Furthermore, the FOXL2 mRNA expression in the mutant group was 4.1-fold that in the wild-type group (Figure 7A), while the FOXL2 protein levels were 2.6-fold that in the wild-type group (Figures 7B,B'). We then assessed the mRNA and protein levels of SRY, FGFR2, FGF9, DMRT1, and FOXL2 genes in the NT2/D1 testicular teratoma cells transfected with the mutant or wild-type MAK3K1 expression plasmid, respectively. The expression of FGFR2 mRNA in the mutant group was 8.5% of that in the wild-type group (Figure 7C), while the protein expression had decreased to 26.5% of that in the wild-type (Figures 7D,D'). The mRNA expression of DMRT1 in the mutant group was 14.1% of that in the wild-type group (Figure 7C), and the protein expression had decreased to 30.5% of the wild-type (Figures 7D,D'). The mRNA and protein expression of FOXL2 gene was 4.1 and 2.4-fold higher, respectively, in the mutant group compared to that in the wild-type group (Figures 7C,D,D'). These results demonstrated that the MAP3K1^R186G variant inhibited the expression of the DMRT1 gene in the testis cells but promoted the expression of the FOXL2 gene in both the testis and ovarian cells. Unexpectedly, the mRNA and protein levels of SRY were 2.6-fold and 2.4-fold higher, respectively, in the mutant group compared to the control group.