Get 20M+ Full-Text Papers For Less Than $1.50/day. Start a 14-Day Trial for You or Your Team.

Learn More →

MYRF haploinsufficiency causes 46,XY and 46,XX disorders of sex development: bioinformatics consideration

MYRF haploinsufficiency causes 46,XY and 46,XX disorders of sex development: bioinformatics... Abstract Disorders of sex development (DSDs) are defined as congenital conditions in which chromosomal, gonadal or anatomical sex is atypical. In many DSD cases, genetic causes remain to be elucidated. Here, we performed a case–control exome sequencing study comparing gene-based burdens of rare damaging variants between 26 DSD cases and 2625 controls. We found exome-wide significant enrichment of rare heterozygous truncating variants in the MYRF gene encoding myelin regulatory factor, a transcription factor essential for oligodendrocyte development. All three variants occurred de novo. We identified an additional 46,XY DSD case of a de novo damaging missense variant in an independent cohort. The clinical symptoms included hypoplasia of Müllerian derivatives and ovaries in 46,XX DSD patients, defective development of Sertoli and Leydig cells in 46,XY DSD patients and congenital diaphragmatic hernia in one 46,XY DSD patient. As all of these cells and tissues are or partly consist of coelomic epithelium (CE)-derived cells (CEDC) and CEDC developed from CE via proliferaiton and migration, MYRF might be related to these processes. Consistent with this hypothesis, single-cell RNA sequencing of foetal gonads revealed high expression of MYRF in CE and CEDC. Reanalysis of public chromatin immunoprecipitation sequencing data for rat Myrf showed that genes regulating proliferation and migration were enriched among putative target genes of Myrf. These results suggested that MYRF is a novel causative gene of 46,XY and 46,XX DSD and MYRF is a transcription factor regulating CD and/or CEDC proliferation and migration, which is essential for development of multiple organs. Introduction Disorders of sex development (DSDs) are defined as congenital conditions in which chromosomal, gonadal or anatomical sex is atypical (1). As per this definition, sex chromosome abnormalities are a type of DSD and include Klinefelter syndrome (47,XXY) with infertility, Turner syndrome (45,XO) with streak gonads, 45,X/46,XY DSD with nearly normal male to Turner syndrome-like phenotype and 46,XX/46,XY DSD with ovotestis. Other types of DSD, namely, 46,XY and 46,XX DSD, are known to be caused by pathogenic variants in specific genes including NR5A1, SRY and SOX9 (2). Recent targeted next-generation sequencing of known mutated genes or whole-exome sequencing (WES) has reported diagnostic yields of 35–40% in DSD cases without sex chromosome aneuploidy (3,4). Thus, in the large proportion of sex chromosome aneuploidy-negative DSD cases, the genetic cause remains unclear. Figure 1 View largeDownload slide Enrichment of rare truncating variants of MYRF in DSD cases. (A, B) Statistical significance (P-value) of gene-based burden of rare damaging variants in DSD cases compared with controls. Two types of damaging variants, LGD model (A) and LGD + D-mis model (B), are shown. Manhattan plot (left) and QQ plot (right) of P-value are shown. Red lines in Manhattan plots, thresholds for statistical significance; red lines in QQ plots, x = y. (C) Four pedigrees with DSD cases harbouring MYRF damaging variants.The properties of the MYRF variants are shown below each individual. Arrows, probands; squares, males; circles, females; closed symbols, affected; open symbols, unaffected; wt, wild-type allele. (D) Poisson distribution of the number of de novo truncating variants of MYRF in 26 samples under the null model of mutation (6). The probability of observing three or more de novo truncating variants of MYRF in 26 individuals (52 alleles) was 1.95E-012. (E) Evolutionary conservation of MYRF Gln443 that is mutated in Patient 5. Note that Gln443 was near an interface for trimerization (11). (F) Hydrogen bonds formed by p.Gln443Pro. The side chain of wild-type Gln443 formed hydrogen bonds with Lys437 and Ala440 (upper), while that of p.Gln443Pro formed two novel hydrogen bonds (pink dotted line) with Ala440. Figure 1 View largeDownload slide Enrichment of rare truncating variants of MYRF in DSD cases. (A, B) Statistical significance (P-value) of gene-based burden of rare damaging variants in DSD cases compared with controls. Two types of damaging variants, LGD model (A) and LGD + D-mis model (B), are shown. Manhattan plot (left) and QQ plot (right) of P-value are shown. Red lines in Manhattan plots, thresholds for statistical significance; red lines in QQ plots, x = y. (C) Four pedigrees with DSD cases harbouring MYRF damaging variants.The properties of the MYRF variants are shown below each individual. Arrows, probands; squares, males; circles, females; closed symbols, affected; open symbols, unaffected; wt, wild-type allele. (D) Poisson distribution of the number of de novo truncating variants of MYRF in 26 samples under the null model of mutation (6). The probability of observing three or more de novo truncating variants of MYRF in 26 individuals (52 alleles) was 1.95E-012. (E) Evolutionary conservation of MYRF Gln443 that is mutated in Patient 5. Note that Gln443 was near an interface for trimerization (11). (F) Hydrogen bonds formed by p.Gln443Pro. The side chain of wild-type Gln443 formed hydrogen bonds with Lys437 and Ala440 (upper), while that of p.Gln443Pro formed two novel hydrogen bonds (pink dotted line) with Ala440. To identify a novel causative gene of DSD, we here performed a case–control WES study comparing the gene-based burden of rare damaging variants between DSD cases (N = 26) and healthy controls (N = 2625). We identified MYRF (MIM: 608329), encoding myelin regulatory factor, as a causative gene of DSD with exome-wide statistical significance. We consequently show convergent evidence indicating a critical role of MYRF in the development of coelomic epithelium (CE)-derived cells (CEDC), which is supported by detailed phenotypic assessment of the DSD cases with the MYRF variants, single-cell gene expression analysis and comprehensive mapping of putative target genes of the transcription factor MYRF. Result In this study, we performed WES on 26 DSD cases and 2625 controls. For the DSD cases, we confirmed that none of the individuals harboured known genetic causes, namely, sex chromosome aneuploidy (Supplementary Material, Fig. S1A and B) or pathogenic/likely pathogenic variants in 68 human-curated genes implicated in DSD (Supplementary Material, Table S1) (2,4). The classification of variant pathogenicity was based on the guidelines from the American College of Medical Genetics and Genomics (5). By using the obtained WES data, we confirmed that there was no apparent population stratification (Supplementary Material, Fig. S1C), no strong relatedness among samples (Supplementary Material, Fig. S1D) and no sample contamination (Supplementary Material, Fig. S1E) in any of the DSD cases and controls, confirming that the quality of our data was sufficient for the following case–control study. We next performed an exome-wide analysis comparing the gene-based burden of rare damaging variants between cases and controls (see Materials and Methods). We analysed two types of damaging variants: likely gene-disrupting variants (LGDs) only or LGD + damaging missense variants (D-mis). For LGD burden, one gene, MYRF at 11q12.2 showed a P-value (8.38E−07) below the threshold for exome-wide statistical significance (2.05E-06, see Materials and Methods) (Fig. 1A, left). For MYRF, there were three LGD variants in cases and none in controls. A quantile–quantile (QQ) plot of observed and expected P-values revealed that there was no inflation in the observed P-values and identified MYRF as the only gene with upward deviation from the theoretical uniform distribution of P-values (Fig. 1A, right). For LGD + D-mis burden, no gene had a P-value below the threshold described above (Fig. 1B, left), and there was no upward outlier from the theoretical distribution of P-value (Fig. 1B, right). Three LGD variants were identified in MYRF (NM_001127392.2): c.278del, p.(Pro93Argfs*7); c.789del, p.(Ser264Alafs*8); and c.1642_1666del, p.(Ala548Thrfs*49) at the fourth, fifth and eighth exons, respectively. All of these variants not located at the last exon or the last 50 bp from the penultimate exon are expected to be subjected to nonsense-mediated mRNA decay, although we could not experimentally test the effect of these LGD variants because cell samples were unavailable in patients. We confirmed the existence of these variants by Sanger sequencing and revealed that all of them are de novo variants (Fig. 1C). In Family 2, the c.789del: p.(Ser264Alafs*8) variant was shared with the affected monozygotic twin. In Family 3, the c.1642_1666del: p.(Ala548Thrfs*49) variant was not found in two unaffected siblings (Fig. 1C). Because all of the LGD variants in MYRF occurred de novo, we additionally investigated whether the observed number of de novo MYRF variants was significantly more than expected in 26 DSD cases under the null model of mutation (6). Under this model, the probability of observing three or more de novo LGD variants in MYRF in 26 individuals was 1.95E-12 (Fig. 1D) and below the threshold of statistical significance (2.05E-06, see Materials and Methods). This result further demonstrates the robust statistical association between LGD variants in MYRF and DSD. The statistical evidence for the enrichment of de novo MYRF variants in DSD motivated us to screen rare MYRF variants in another DSD cohort. Of 68 DSD cases of unknown genetic cause, we identified a de novo rare MYRF missense variant, c.1328A>C: p.(Gln443Pro), in Family 4 (Fig. 1C). The pathogenicity of this variant was supported by multiple lines of evidence as follows. (1) Most in silico tools predicted the variant to be damaging (Supplementary Material, Table S2) (7–10). (2) The amino acid at the variant was conserved among species (Fig. 1E). (3) There is no missense variant affecting the same amino acid (Gln443) in general populations: ESP6500, gnomAD and ToMMo3.5KJPN. (4) The variant was located adjacent to the interface for homo-trimerization (11) (Fig. 1E), where missense variants are depleted in general populations, as shown in the DECIPHER browser (12) (data not shown). The missense variant was predicted by Swiss PDB viewer structural analysis (Protein Data Bank ID: 5H5P) to intramolecularly form a novel hydrogen bond (Fig. 1F) and might change the structure of the adjacent interface for homo-trimerization. When we analysed the probability of observing four or more de novo LGD or missense variants in MYRF among 94 (26 + 68) individuals under the null model of mutation (6), we calculated the probability as 1.0E-10, which was again below the threshold of exome-wide statistical significance (2.05E-06). Table 1 Summary of clinical and genetic findings in DSD cases with MYRF damaging variants Case Pt1 Monozygotic twin Pt4 Pt5 Pt2 Pt3  Karyotype 46,XY 46,XX 46,XX 46,XY 46,XY  Social sex F—>M at 4 y F F M M  Present age 14 years 21 years 21 years 4 years 4 months 2 years 5 months  Genital findings   Age at examination 2 months 15 years 15 years 3 momths 0 (At birth)   Tanner stage I III (breast); II (pubis) II (breast); II (pubis) Prepubertal Prepubertal   Testis size (mL) Unknown - - Small (<1 mL, L) 0.9 × 0.5 × 0.6 cm on echo Small (<1 mL, L and R)   Cryptorchidism Abdominal (R); inguinal (L) - - Abdominal (L and R) Inguinal (L and R)   Vanishing testis No - - + (R) No   Hypospadias Perineal type - - Penile type Penoscrotal type   Microphallus 2.3 cm (−3.0 SD) - - 1.8 cm (−3.0 SD) 1.5 cm (<−3.0 SD)   Abnormal scrotum Anterior-positioned and bifid - - No Small and bifid   Ovary - Absent (L); small (R) Absent (L); small (R) - -   Fallopian tube - Absent Absent - -   Uterus - Absent (L); restiform cervix (R) Absent (L); restiform cervix (R) - -   Vagina - Absent, surrounded by a cyst Absent, surrounded by a cyst - -  Endocrine findings   Age at examination 2 months 15 years 15 years 1 year 0 month 3 months 7 months   LH (IU/L) baseline NA 45.2 25.7 0.73 (0.2–1.9) 8.9 (0.2–1.9) 0.3 (0.2–1.9)   LH (IU/L) GnRH-stimulated NA NA NA 9.49 (1.1–6.0) NA 8.4 (1.1–6.0)   FSH (IU/L) baseline NA 70.1 113.1 3.91 (<0.3–2.4) 8.9 (<0.3–2.4) 3.4 (<0.3–2.4)   FSH (IU/L) GnRH-stimulated NA NA NA 21.62 (1.9–7.6) NA 23.9 (1.9–7.6)   T (ng/dL) baseline 125 (120 ~ 400) 55.0 54.9 <5 (<5) 116 (120–400) <5 (<5)   T (ng/dL) hCG-stimulated 440 (> 200) NA NA 518 (> 200) NA 78 (>200)   DHT (ng/dL) baseline NA NA NA <5 (<5) 20 NA   AMH (ng/mL) baseline NA NA NA 18.3 (43.3–79.3) 19.6 (43.3–79.3) 15.0 (43.3–79.3)  Treatment   Orchidopexy (age) 1 years - - 1 year 7 months Not yet   Genitoplasty (age) 2 years 15 years 15 years 2 years 11 month Not yet   TE (25 mg i.m.) (age) 3 × (6 ~ 10 months) - - No 1 × (5 months)   DHT cream (age) No - - Yes (2 years 8 months) No   PL increment (cm) 2.3→3.9 (+ 1.1 SD) - - Increased (no record) 1.5→3.0 (+1.2 SD)  Central nervous system   Developmental delay No No No Yes No   DQ NA NA NA 63 NA   Seizure No No No No No   ASD/ADHD No No No NA No  Others   Diaphragmatic hernia No No No Yes (L) No Case Pt1 Monozygotic twin Pt4 Pt5 Pt2 Pt3  Karyotype 46,XY 46,XX 46,XX 46,XY 46,XY  Social sex F—>M at 4 y F F M M  Present age 14 years 21 years 21 years 4 years 4 months 2 years 5 months  Genital findings   Age at examination 2 months 15 years 15 years 3 momths 0 (At birth)   Tanner stage I III (breast); II (pubis) II (breast); II (pubis) Prepubertal Prepubertal   Testis size (mL) Unknown - - Small (<1 mL, L) 0.9 × 0.5 × 0.6 cm on echo Small (<1 mL, L and R)   Cryptorchidism Abdominal (R); inguinal (L) - - Abdominal (L and R) Inguinal (L and R)   Vanishing testis No - - + (R) No   Hypospadias Perineal type - - Penile type Penoscrotal type   Microphallus 2.3 cm (−3.0 SD) - - 1.8 cm (−3.0 SD) 1.5 cm (<−3.0 SD)   Abnormal scrotum Anterior-positioned and bifid - - No Small and bifid   Ovary - Absent (L); small (R) Absent (L); small (R) - -   Fallopian tube - Absent Absent - -   Uterus - Absent (L); restiform cervix (R) Absent (L); restiform cervix (R) - -   Vagina - Absent, surrounded by a cyst Absent, surrounded by a cyst - -  Endocrine findings   Age at examination 2 months 15 years 15 years 1 year 0 month 3 months 7 months   LH (IU/L) baseline NA 45.2 25.7 0.73 (0.2–1.9) 8.9 (0.2–1.9) 0.3 (0.2–1.9)   LH (IU/L) GnRH-stimulated NA NA NA 9.49 (1.1–6.0) NA 8.4 (1.1–6.0)   FSH (IU/L) baseline NA 70.1 113.1 3.91 (<0.3–2.4) 8.9 (<0.3–2.4) 3.4 (<0.3–2.4)   FSH (IU/L) GnRH-stimulated NA NA NA 21.62 (1.9–7.6) NA 23.9 (1.9–7.6)   T (ng/dL) baseline 125 (120 ~ 400) 55.0 54.9 <5 (<5) 116 (120–400) <5 (<5)   T (ng/dL) hCG-stimulated 440 (> 200) NA NA 518 (> 200) NA 78 (>200)   DHT (ng/dL) baseline NA NA NA <5 (<5) 20 NA   AMH (ng/mL) baseline NA NA NA 18.3 (43.3–79.3) 19.6 (43.3–79.3) 15.0 (43.3–79.3)  Treatment   Orchidopexy (age) 1 years - - 1 year 7 months Not yet   Genitoplasty (age) 2 years 15 years 15 years 2 years 11 month Not yet   TE (25 mg i.m.) (age) 3 × (6 ~ 10 months) - - No 1 × (5 months)   DHT cream (age) No - - Yes (2 years 8 months) No   PL increment (cm) 2.3→3.9 (+ 1.1 SD) - - Increased (no record) 1.5→3.0 (+1.2 SD)  Central nervous system   Developmental delay No No No Yes No   DQ NA NA NA 63 NA   Seizure No No No No No   ASD/ADHD No No No NA No  Others   Diaphragmatic hernia No No No Yes (L) No NOTE: LH indicates luteinizing hormone; i.m., intramuscular injection; FSH, follicle-stimulating hormone; DHT, dihydrotestosterone; TE, testosterone enanthate; PL, penile length; DQ, development quotient measured with Kyoto Scale of Psychological Development; ASD, autism spectrum disorder; ADHD, attention-deficit hyperactivity disorder. View Large Table 1 Summary of clinical and genetic findings in DSD cases with MYRF damaging variants Case Pt1 Monozygotic twin Pt4 Pt5 Pt2 Pt3  Karyotype 46,XY 46,XX 46,XX 46,XY 46,XY  Social sex F—>M at 4 y F F M M  Present age 14 years 21 years 21 years 4 years 4 months 2 years 5 months  Genital findings   Age at examination 2 months 15 years 15 years 3 momths 0 (At birth)   Tanner stage I III (breast); II (pubis) II (breast); II (pubis) Prepubertal Prepubertal   Testis size (mL) Unknown - - Small (<1 mL, L) 0.9 × 0.5 × 0.6 cm on echo Small (<1 mL, L and R)   Cryptorchidism Abdominal (R); inguinal (L) - - Abdominal (L and R) Inguinal (L and R)   Vanishing testis No - - + (R) No   Hypospadias Perineal type - - Penile type Penoscrotal type   Microphallus 2.3 cm (−3.0 SD) - - 1.8 cm (−3.0 SD) 1.5 cm (<−3.0 SD)   Abnormal scrotum Anterior-positioned and bifid - - No Small and bifid   Ovary - Absent (L); small (R) Absent (L); small (R) - -   Fallopian tube - Absent Absent - -   Uterus - Absent (L); restiform cervix (R) Absent (L); restiform cervix (R) - -   Vagina - Absent, surrounded by a cyst Absent, surrounded by a cyst - -  Endocrine findings   Age at examination 2 months 15 years 15 years 1 year 0 month 3 months 7 months   LH (IU/L) baseline NA 45.2 25.7 0.73 (0.2–1.9) 8.9 (0.2–1.9) 0.3 (0.2–1.9)   LH (IU/L) GnRH-stimulated NA NA NA 9.49 (1.1–6.0) NA 8.4 (1.1–6.0)   FSH (IU/L) baseline NA 70.1 113.1 3.91 (<0.3–2.4) 8.9 (<0.3–2.4) 3.4 (<0.3–2.4)   FSH (IU/L) GnRH-stimulated NA NA NA 21.62 (1.9–7.6) NA 23.9 (1.9–7.6)   T (ng/dL) baseline 125 (120 ~ 400) 55.0 54.9 <5 (<5) 116 (120–400) <5 (<5)   T (ng/dL) hCG-stimulated 440 (> 200) NA NA 518 (> 200) NA 78 (>200)   DHT (ng/dL) baseline NA NA NA <5 (<5) 20 NA   AMH (ng/mL) baseline NA NA NA 18.3 (43.3–79.3) 19.6 (43.3–79.3) 15.0 (43.3–79.3)  Treatment   Orchidopexy (age) 1 years - - 1 year 7 months Not yet   Genitoplasty (age) 2 years 15 years 15 years 2 years 11 month Not yet   TE (25 mg i.m.) (age) 3 × (6 ~ 10 months) - - No 1 × (5 months)   DHT cream (age) No - - Yes (2 years 8 months) No   PL increment (cm) 2.3→3.9 (+ 1.1 SD) - - Increased (no record) 1.5→3.0 (+1.2 SD)  Central nervous system   Developmental delay No No No Yes No   DQ NA NA NA 63 NA   Seizure No No No No No   ASD/ADHD No No No NA No  Others   Diaphragmatic hernia No No No Yes (L) No Case Pt1 Monozygotic twin Pt4 Pt5 Pt2 Pt3  Karyotype 46,XY 46,XX 46,XX 46,XY 46,XY  Social sex F—>M at 4 y F F M M  Present age 14 years 21 years 21 years 4 years 4 months 2 years 5 months  Genital findings   Age at examination 2 months 15 years 15 years 3 momths 0 (At birth)   Tanner stage I III (breast); II (pubis) II (breast); II (pubis) Prepubertal Prepubertal   Testis size (mL) Unknown - - Small (<1 mL, L) 0.9 × 0.5 × 0.6 cm on echo Small (<1 mL, L and R)   Cryptorchidism Abdominal (R); inguinal (L) - - Abdominal (L and R) Inguinal (L and R)   Vanishing testis No - - + (R) No   Hypospadias Perineal type - - Penile type Penoscrotal type   Microphallus 2.3 cm (−3.0 SD) - - 1.8 cm (−3.0 SD) 1.5 cm (<−3.0 SD)   Abnormal scrotum Anterior-positioned and bifid - - No Small and bifid   Ovary - Absent (L); small (R) Absent (L); small (R) - -   Fallopian tube - Absent Absent - -   Uterus - Absent (L); restiform cervix (R) Absent (L); restiform cervix (R) - -   Vagina - Absent, surrounded by a cyst Absent, surrounded by a cyst - -  Endocrine findings   Age at examination 2 months 15 years 15 years 1 year 0 month 3 months 7 months   LH (IU/L) baseline NA 45.2 25.7 0.73 (0.2–1.9) 8.9 (0.2–1.9) 0.3 (0.2–1.9)   LH (IU/L) GnRH-stimulated NA NA NA 9.49 (1.1–6.0) NA 8.4 (1.1–6.0)   FSH (IU/L) baseline NA 70.1 113.1 3.91 (<0.3–2.4) 8.9 (<0.3–2.4) 3.4 (<0.3–2.4)   FSH (IU/L) GnRH-stimulated NA NA NA 21.62 (1.9–7.6) NA 23.9 (1.9–7.6)   T (ng/dL) baseline 125 (120 ~ 400) 55.0 54.9 <5 (<5) 116 (120–400) <5 (<5)   T (ng/dL) hCG-stimulated 440 (> 200) NA NA 518 (> 200) NA 78 (>200)   DHT (ng/dL) baseline NA NA NA <5 (<5) 20 NA   AMH (ng/mL) baseline NA NA NA 18.3 (43.3–79.3) 19.6 (43.3–79.3) 15.0 (43.3–79.3)  Treatment   Orchidopexy (age) 1 years - - 1 year 7 months Not yet   Genitoplasty (age) 2 years 15 years 15 years 2 years 11 month Not yet   TE (25 mg i.m.) (age) 3 × (6 ~ 10 months) - - No 1 × (5 months)   DHT cream (age) No - - Yes (2 years 8 months) No   PL increment (cm) 2.3→3.9 (+ 1.1 SD) - - Increased (no record) 1.5→3.0 (+1.2 SD)  Central nervous system   Developmental delay No No No Yes No   DQ NA NA NA 63 NA   Seizure No No No No No   ASD/ADHD No No No NA No  Others   Diaphragmatic hernia No No No Yes (L) No NOTE: LH indicates luteinizing hormone; i.m., intramuscular injection; FSH, follicle-stimulating hormone; DHT, dihydrotestosterone; TE, testosterone enanthate; PL, penile length; DQ, development quotient measured with Kyoto Scale of Psychological Development; ASD, autism spectrum disorder; ADHD, attention-deficit hyperactivity disorder. View Large We next investigated the phenotype of the DSD cases with de novo MYRF variants in the four families (Table 1). The three XY patients, Patients 1, 4 and 5, had micro-penis, hypospadias, small testis and cryptorchidism (Supplementary Material, Fig. S2). Laboratory data showed low levels of testosterone (T) and anti-Müllerian hormone (AMH), as well as hypergonadotropism. Intramuscular injection of testosterone enanthate effectively increased the penile length of Patient 4. The low level of AMH is indicative of defects in Sertoli cells. The other findings can be caused by a low level of T during the foetal period, which is cooperatively produced by Sertoli and Leydig cells. One patient showed congenital diaphragmatic hernia (CDH), a cause of which is a defect in the development of pleuroperitoneal folds (13). In the monozygotic twin XX patients of Family 2, uterus, fallopian tube and ovary were hypoplastic, and laboratory data showed hypergonadotropism. These phenotypes can be explained by the defective development of Müllerian derivatives and ovaries. Therefore, the majority of the symptoms observed in the DSD cases of MYRF pathogenic variants can be attributable to defects in the following cells and tissues: Müllerian derivatives, ovaries, Sertoli cells, Leydig cells and pleuroperitoneal folds. Interestingly, all of these cells/tissues are CEDC or partly consist of CEDC (Fig. 2A) and developed via processes of proliferaiton, migration and epithelial-to-mesenchymal transition (EMT) and mesenchymal-to-epithelial transition (MET) (14–18). Thus, the phenotype of the patients suggested that MYRF might play a critical role in the development of CEDC, such as their proliferation, migration, EMT and/or MET. We next investigated the expression patterns of MYRF in foetal gonads. To this end, we analysed public single-cell RNA sequencing (scRNA-seq) data of the human male and female aorta–gonad–mesonephros (AGM) region (4 to 5 wpf) and foetal gonads (7 to 26 wpf) (19). By analysing female gonad data by t-Distributed Stochastic Neighbour Embedding (t-SNE), we identified 11 clusters of possibly different cell identities (C1–11, Fig. 2B). We found that the expression of MYRF was high in C9, 10 and 11 (Fig. 2C and D). We determined the cell identities of each cluster mainly focusing on C9, 10 and 11 as follows (Supplementary Material, Fig. 3). Some cells of C9, C10 and C11 expressed UPK3B, a CE marker (Fig. 2D) (20). C9, C10 and C11 cells also expressed FOXL2, a medullary granulosa cell lineage marker (21), and NR2F2, a putative thecal cell lineage marker (22) (Fig. 2C and D). C9, C10 and C11 rarely expressed LGR5, a cortical granulosa cell lineage marker (22), possibly due to an insufficient sequencing depth (data not shown). Information about CCDC182, another granulosa cell lineage marker (23), was not included in the data, possibly due to the use of old gene annotations. Thus, C9, C10 and C11 could be mixtures of CE, cells of the granulosa cell lineage and cells of the thecal cell lineage. Upon the temporal resolution of each cluster, C9, C10 and C11 mainly consisted of cells at 5–10, 14 and 18–26 wpf and were considered as early, middle and late cells, respectively (Supplementary Material, Fig. S3A). A few cells in the upper part of C9 showed PECAM1 (24) and CDH5 (25) positivity and were considered as endothelial cells (Supplementary Material, Fig. S3M and N). Taking these results together, the properties of cells with high MYRF expression were compatible with CE and CEDC in females, that is, cells of granulosa cell and putative theca cell lineage (Fig. 2A). Figure 2 View largeDownload slide High expression of MYRF in CE, granulosa cell lineage and/or putative theca cell lineage in female gonads. (A) Schematic representation of the relationships among CE and CEDC in foetal gonads. The origin of Leydig cell lineage is partly CE while that of theca cell lineage is not known. Previously reported gene markers for each cell type are as follows: UPK3B (20), FOXL2 (21), NR2F2 (22,27), LGR5 (22), CCDC182 (23), WT1 (26), AMH (26), ARX (28) and CYP17A1 (55). (B–D) scRNA-seq analysis of human female foetal gonads. (B) t-SNE plots showing identities of each cell cluster. Each dot indicates each cell. Dots are coloured based on cell identities determined by t-SNE analysis. Cell type of each cluster was determined by gene markers (C and D) and wpf (Supplementary Material, Fig. S3A). (C) t-SNE plots showing log-normalized expression level of MYRF and marker genes as shown on the right side. Analysed gene is shown above each plot. (D) Bar plot of log-normalized expression level of MYRF and marker genes as shown above each plot. Each bar indicates each cell. Clusters are shown in different colours, so that each cluster is differentiated in the graph. Figure 2 View largeDownload slide High expression of MYRF in CE, granulosa cell lineage and/or putative theca cell lineage in female gonads. (A) Schematic representation of the relationships among CE and CEDC in foetal gonads. The origin of Leydig cell lineage is partly CE while that of theca cell lineage is not known. Previously reported gene markers for each cell type are as follows: UPK3B (20), FOXL2 (21), NR2F2 (22,27), LGR5 (22), CCDC182 (23), WT1 (26), AMH (26), ARX (28) and CYP17A1 (55). (B–D) scRNA-seq analysis of human female foetal gonads. (B) t-SNE plots showing identities of each cell cluster. Each dot indicates each cell. Dots are coloured based on cell identities determined by t-SNE analysis. Cell type of each cluster was determined by gene markers (C and D) and wpf (Supplementary Material, Fig. S3A). (C) t-SNE plots showing log-normalized expression level of MYRF and marker genes as shown on the right side. Analysed gene is shown above each plot. (D) Bar plot of log-normalized expression level of MYRF and marker genes as shown above each plot. Each bar indicates each cell. Clusters are shown in different colours, so that each cluster is differentiated in the graph. Figure 3 View largeDownload slide High expression of MYRF in Sertoli and Leydig cell lineages in male gonads. scRNA-seq analysis of human male foetal gonads. (A) t-SNE plots showing identities of each cell cluster. Each dot indicates each cell. Dots are coloured based on cell identities determined by t-SNE analysis. Cell type of each cluster was determined by gene markers (B, C) and wpf (Supplementary Material, Fig. S4A). (B) t-SNE plots showing log-normalized expression level of MYRF and marker genes as shown on the right side. Analysed gene is shown above each splot. (C) Bar plot of log-normalized expression level of MYRF and marker genes as shown above each plot. Each bar indicates each cell. Clusters are shown in different colours so that each cluster is differentiated in the graph. In addition, early (red) and late (yellow green) Sertoli cells in C6; early (green), middle (blue) and late (purple) Leydig cells in C7; and migrating (red) and other FGC (yellow green) in C1 are coloured differently within each cluster. Figure 3 View largeDownload slide High expression of MYRF in Sertoli and Leydig cell lineages in male gonads. scRNA-seq analysis of human male foetal gonads. (A) t-SNE plots showing identities of each cell cluster. Each dot indicates each cell. Dots are coloured based on cell identities determined by t-SNE analysis. Cell type of each cluster was determined by gene markers (B, C) and wpf (Supplementary Material, Fig. S4A). (B) t-SNE plots showing log-normalized expression level of MYRF and marker genes as shown on the right side. Analysed gene is shown above each splot. (C) Bar plot of log-normalized expression level of MYRF and marker genes as shown above each plot. Each bar indicates each cell. Clusters are shown in different colours so that each cluster is differentiated in the graph. In addition, early (red) and late (yellow green) Sertoli cells in C6; early (green), middle (blue) and late (purple) Leydig cells in C7; and migrating (red) and other FGC (yellow green) in C1 are coloured differently within each cluster. We repeated the above analysis for public data of scRNA-seq on human male foetal gonad. We performed t-SNE analysis and identified eight clusters of possibly different cell identities (C1–8, Fig. 3A). We found that the expression of MYRF was high in C6 diffusely and in some cells of C7 (Fig. 3B). We determined the cell identities of each cluster, mainly focusing on C6 and C7, as follows (Fig. 3B and C and Supplementary Material, Fig. S4) (19). All clusters expressed UPK3B similarly and at a low level (data not shown), possibly indicating that the data contained no CE. C6 expressed WT1 (26) and AMH (26) corresponding to the Sertoli cell lineage (Fig. 3B and C). C7 expressed NR2F2 (27) and ARX (28) and corresponded to the Leydig cell lineage (Fig. 3B and C). Because a sub-cluster of the right part of C6 was mainly from before 10 wpf (Supplementary Material, Fig. S4A), this sub-cluster was considered to consist early cells of the Sertoli cell lineage, while the other part was probably late cells of this lineage. Because a sub-cluster of the upper part of C7 was mainly from before 12 wpf (Supplementary Material, Fig. S4A), the sub-cluster was considered as early cells of the Leydig cell lineage. The few cells at the left part showed positivity for CYP17A1 (29) (Fig. 3B and C), a marker for late cells of the Leydig cell lineage. Thus, the properties of cells with high MYRF expression were compatible with CEDC in males, namely, cells of the Sertoli and Leydig cell lineages (Fig. 2A). Figure 4 View largeDownload slide GO enrichment analysis of putative target genes of Myrf. GO enrichment analysis of putative Myrf target genes in rat differentiating oligodendrocyte progenitor cells. (A) Bar plot of statistical significance (PBH) for enrichment of GO terms of interest related to proliferation, migration and EMT/MET. Bars are coloured according to the categories of GO terms as follows: white for proliferation, grey for migration and black for EMT/MET. The dotted horizontal line indicates PBH of 0.05. (B) Functions of clusters of GO terms with a raw P-value of <0.05 with overlapping genes. Node (white to red circle), GO; node size, correlated with the number of contained genes; node colour, correlated with raw P-value shown as a colour bar at the bottom left; edge (blue lines), overlap of contained genes between GO terms; black oval, GO terms that were similar to each other. Figure 4 View largeDownload slide GO enrichment analysis of putative target genes of Myrf. GO enrichment analysis of putative Myrf target genes in rat differentiating oligodendrocyte progenitor cells. (A) Bar plot of statistical significance (PBH) for enrichment of GO terms of interest related to proliferation, migration and EMT/MET. Bars are coloured according to the categories of GO terms as follows: white for proliferation, grey for migration and black for EMT/MET. The dotted horizontal line indicates PBH of 0.05. (B) Functions of clusters of GO terms with a raw P-value of <0.05 with overlapping genes. Node (white to red circle), GO; node size, correlated with the number of contained genes; node colour, correlated with raw P-value shown as a colour bar at the bottom left; edge (blue lines), overlap of contained genes between GO terms; black oval, GO terms that were similar to each other. We also analysed Myrf expression in another dataset, the Mouse Cell Atlas (30). This dataset contains scRNA-seq data of male and female mouse gonads at E14.5, and we analysed them as we did for the above human foetal gonad data (Supplementary Material, Fig. S5 and S6). In the male and female mouse gonad data, Myrf was rarely expressed, possibly due to a low sequencing depth (Supplementary Material, Fig. S5B and S6B). The cluster in which the mean Myrf expression was the highest was CE (C9 in Supplementary Material, Figs S5 and S6, C8) in both male and female foetal gonads (Supplementary Material, Figs S5C, D and 6C, D). Other clusters expressing Myrf were cells of the granulosa cell lineage and putative thecal cell lineage and FGC in female foetal gonad, while they were early cells of the Leydig cell lineage and cells of the Sertoli cell lineage in male foetal gonad (Supplementary Material, Figs. S5C, D and 6C, D). Myelin regulatory factor protein encoded by MYRF is known as a transcription factor positively regulating gene expression (31,32). We therefore analysed putative target genes of MYRF by using public data of ChIP-seq for Myrf in rat differentiating oligodendrocyte progenitor cells (31). The experiment was performed in a type of cell that differs from those in which we are interested. However, because target genes of transcription factors largely overlap between different types of cell (33), we could partially reveal the target genes of MYRF in affected cells in our patients. We found that Myrf peaks identified by ChIP-seq are highly condensed in regions proximal to transcription start site (TSS; highest within ±10 kb; Supplementary Material, Fig. S7). The results were consistent with a previously reported finding (31) and suggest that these proximal genes are putative targets of Myrf. Thus, we performed a Gene Ontology (GO) enrichment analysis of putative Myrf-target genes whose TSS is close to Myrf peaks (<10 kb) using CHIP-ENRICH (number of putative target genes = 478). As described above, the available ChIP-seq data are derived from oligodendrocyte progenitor cells, and it is expected the target genes are enriched for oligodendritic and other neural processes. Therefore, we first selected a total of 31 GO terms of interest related to proliferation, migration and EMT/MET prior to the analysis (see Materials and Methods) and performed a CHIP-ENRICH analysis. We found that eight and one GO terms for migration and proliferation were significantly [Benjamini–Hochberg-adjusted P-value (PBH) < 0.05] enriched among the putative Myrf target genes, respectively (Fig. 4A). Next, we performed an exploratory enrichment analysis of all GO terms. We identified only three terms with PBH < 0.05: negative regulation of cellular component organization (GO: 0051129), negative regulation of microtubule polymerization or depolymerization (GO: 0031111) and regulation of receptor binding (GO: 1900120) (Supplementary Material, Table S3). To obtain broader insights into the functions of putative Myrf-target genes, we applied a more permissive threshold, raw P < 0.05, and obtained 309 GO terms with nominal significance (Supplementary Material, Table S3). To summarize the functions of these GO terms, we grouped the terms into clusters reflecting their overlapping genes. The two largest clusters were, as expected, related to brain development and development of the nervous system, followed by the regulation of cell migration, regulation of cell-matrix adhesion, regulation of microtubule organization, regulation of response to reactive oxygen species (especially cell death), regulation of JAK-STAT cascade, regulation of kinase (especially MAPK cascade), development of blood vessels and embryonic development (Fig. 4B). Discussion In this study, we identified MYRF as a strong candidate for a gene causative of DSD, for which robust statistically significant results were obtained. The pathogenicity of MYRF damaging variants was supported by multiple lines of evidence: (1) exome-wide significant enrichment of rare LGD variants at MRYF in DSD cases, (2) exome-wide significant enrichment of de novo truncating variants at MYRF in DSD cases, (3) high and selective expression of MYRF in cells and tissues affected in the DSD cases with MYRF pathogenic variants and (4) enrichment of genes related to proliferation and migration among putative target genes of Myrf. Moreover, Pinz et al. (34) recently reported two male cases of DSD and CDH with MYRF truncating variants. While their study alone could not provide robust statistical evidence due to the small number of identified mutation carriers and the lack of a comparison with a control dataset, together with our exome-wide significant findings, there is convincing evidence that MYRF haploinsufficiency causes a type of DSD. Our patients with MYRF variants showed not only gonadal hypoplasia in both 46,XY and 46,XX but also Müllerian duct hypoplasia in 46,XX subjects. The developments of indifferent gonad and Müllerian duct are different biological processes, and gonadal agenesis usually permits Müllerian development. However, agenesis of both gonad and Müllerian duct has been reported in patients with 46,XY (35), 46,XX (36) and 46,X,i(Xq) (37) (since such 46,XY patients with agenesis of both gonad and Müllerian duct exhibit complete female genitalia, it is unlikely that testis tissue is once formed to produce AMH and then regresses towards gonadal hypoplasia). These cases in our study could support the presence of a gene involved in development of both indifferent gonad and Müllerian duct, and the current cases in our study suggest that MYRF is indeed such a candidate. According to the results of our bioinformatic analyses utilizing scRNA-seq and Myrf ChIP-seq data, we found that MYRF/Myrf are highly expressed in CE/CEDC in both human and mouse foetal gonads. GO enrichment analysis of putative Myrf target genes revealed enrichment of GO terms related to migration and proliferation among them. Interpreting these results in an integrated manner, it can be speculated that loss of function of a copy of the MYRF transcription factor leads to the transcriptional dysregulation of genes related to migration and proliferation in CE/CEDC, and then causes the observed phenotypes that can be mostly explained by defects in CE/CEDC development. Nevertheless, this hypothesis could be too simple, and further details of pathomechanisms should be investigated by other approaches, including an analysis of living animals such as heterozygous Myrf knockout mice. We also need to take into account that we used Myrf ChIP-seq data derived from oligodendrocytes (31), which is not a cell type directly affected in our DSD cases. While in this study we report strong evidence indicating a causal relationship between MYRF haploinsufficiency and DSD, MYRF was originally identified as a transcriptional regulator essential for the development and maintenance of central nervous system myelination (32,38). Indeed, MYRF has been reported to be a causative gene of encephalitis/encephalopathy with a reversible splenial lesion (MERS) and MERS-like diseases triggered by febrile illness (39). Considering the expression levels of MYRF in oligodendrocytes and foetal gonads, we found that the Myrf expression level [log10 (mean TPM in cluster +1)] in epithelial cell clusters, which correspond to CE clusters in our analyses, was 4.7 in male foetal gonad and 3.9 in female foetal gonad (data of Mouse Cell Atlas), which was comparable to that in the oligodendrocyte precursor cell cluster (4.8, data of Mouse Cell Atlas). Thus, MYRF can play an important role in CE/CEDC, and it would not be surprising that MYRF damaging variants cause abnormal sex development. Regarding the phenotype–genotype relationship, three out of four DSD cases in our study, and two of two patients in the study of Pinz et al. (34), carry truncating de novo mutations in MYRF. In contrast, the variant reported in families with MERS-like diseases was a missense variant (p.Gln403Arg) that partially decreases the transcriptional activity in vitro (39). Therefore, partial disruption of MYRF function may cause MERS-like diseases that only become apparent after febrile illness and full haploinsufficiency of MYRF would cause DSD. In addition, while in our DSD cases abnormalities of the nervous system were not observed, except in Patient 5 with developmental delay, symptoms that would be induced by febrile illness should be carefully followed. In summary, in this study we demonstrate a causal relationship between MYRF truncating variants and DSD, supported by robust statistical evidence. Our findings can contribute to an improved understanding and diagnosis of DSD. In addition, our integrative bioinformatic analyses provide insights into the putative mechanism by which MYRF haploinsufficiency leads to the observed DSD phenotypes. Finally, the results of this study expand the spectrum of diseases caused by MYRF mutations and clearly point out that MYRF plays important roles not only in the differentiation and maintenance of oligodendrocytes but also in sex development. Materials and Methods Study participants and ethics Informed consent was obtained from all of the study participants. When participants were younger than 20 years, we obtained informed consent from their parents. This study was approved by the institutional review board of Yokohama City University School of Medicine (Yokohama, Japan). Trained obstetricians, gynaecologists and/or paediatricians performed the DSD diagnosis and the detailed clinical assessment of the patients. Control samples in the following case–control analysis were from unrelated healthy parents who had children suffering from rare diseases. Whole-exome sequencing WES was performed for 26 DSD and 2625 control samples as previously described (37–40). In brief, genomic DNA was captured with a SureSelect Human All Exon V4, V5 or V6 Kit (Agilent Technologies, Santa Clara, CA, USA): V4 in 8, V5 in 16 and V6 in 2 in 26 DSD samples while V4 in 144, V5 in 1567 and V6 in 914 samples of 2625 control samples. Sequencing was performed on an Illumina HiSeq 2500 (Illumina, San Diego, CA, USA) with 101 bp paired-end reads. Reads were aligned to the human reference genome (GRCh37.1/hg19) using Novoalign v3.02.13. Polymerase chain reaction (PCR) duplicates were removed using Picard. Local realignments around indels and base quality score recalibration were performed with the Genome Analysis Toolkit (GATK) 3.7-0, and analysis-ready binary alignment map (BAM) files were generated (41). From the BAM files, variants in all case and control samples were jointly genotyped with GATK HaplotypeCaller and hard-filtered in accordance with GATK Best Practice (42). To remove possible bias from different proportions of capture kits used in DSD and control groups, we selected variants only at bait regions overlapped among V4, V5 and V6 capture kits with a margin of 100 bp using merge and intersect function of bedtools. Variants were annotated using dbNSFP in ANNOVAR and snpEff (43–45). Case–control analysis of rare variant burden The burden of rare damaging variants at each protein-coding gene was compared between 26 DSD cases and 2625 controls. We defined damaging variants in the following two ways: LGD or LGD + D-mis at each gene. We defined variants of HIGH impact in snpEff (e.g. stop_gained, frameshift_variant and splice_donor/acceptor_variant) as LGD and missense variants with a SIFT (8) score of <0.05 and a PolyPhen2 HVAR (7) score of >0.90 as D-mis. We defined rare variants as follows: (1) minor allele frequency (MAF) of <0.0005 in databases of variants in general populations: ESP6500, gnomAD (46) and ToMMo3.5KJPN (47) and (2) not existing in our intramural databases of 575 Japanese controls. We compared the burden of these variants in each gene between DSD cases and controls by a collapsing method (48) followed by Fisher’s exact test using RVTESTS (49). We conservatively defined the threshold P-value for statistical significance by applying Bonferroni correction with the total number of tests in our LGD only and LGD + D-mis models. After excluding genes with no rare damaging variant among the overall cohort (26 DSD cases and 2625 controls), to which Fisher’s exact test cannot be applied, there were 9217 and 15 177 testable genes in LGD only and LGD + D-mis models, respectively. With these numbers of testable genes, we defined the threshold for exome-wide significance as 2.05E−06 [= 0.05/(9217 + 15 177)]. Analysis of burden of de novo variants The burden of the de novo truncating variant at MYRF in 26 DSD cases was analysed under a null model of mutation (6) that compares observed and expected numbers of de novo variants in a gene of interest. Briefly, numbers of de novo variants are expected to follow a Poisson distribution with λ (expected number of de novo variants) calculated as follows: (rate of de novo truncating variants at gene of interest per allele per generation) × (number of tested alleles). Samocha et al. (6) reported the probabilities of de novo nonsense, frameshift, canonical splice-site and missense variants per generation in a copy of MYRF as 1.12E-06, 2.19E-06, 1.06E-06 and 3.29E-05 and thereby the probability of LGD variants or LGD + missense variants can be considered as the sum of these probabilities (4.36E-06 or 3.73E-05, respectively). Thus, the expected count of de novo MYRF LGD variants in 26 cases or LGD + missense variants in 94 cases is expected to follow a Poisson distribution with a mean of 2.27E-04 or 7.01E-03, respectively (= 4.36E-06 × 26 × 2 or 3.73E-05 × 94 × 2). Because we observed three de novo LGD variants in MYRF in 26 DSD cases and four de novo LGD + missense variants in 94 DSD cases, the probability of three or more counts in the Poisson distribution was calculated with the POISSON function of Microsoft Excel. The threshold for statistical significance was 2.05E-06 as described in the above case–control analysis of the rare variant burden. Sanger sequencing The identified MYRF variants in the probands and the family members were analysed by Sanger sequencing using standard methods. Briefly, PCR products were purified with alkaline phosphatase and exonuclease 1 and sequenced with BigDye Terminator v3.1 Cycle Sequencing kit (Applied Biosystems, Foster City, CA, USA) on a 3500 DNA Sequencing Analyzer (Life Technologies, Carlsbad, CA, USA). Primer sequences are available on request. scRNA-seq analysis We reanalysed two public sets of data of scRNA-seq: (1) data of the human AGM region at 5 weeks post-fertilization (wpf) and foetal gonads from 7 to 26 wpf (19) and (2) data of mouse foetal gonads (30). We analysed the data with Seurat (50) in accordance with the developer’s guide. Briefly, we retained genes expressed in at least three cells and retained cells with 200–2500 detected genes using the FilterCells function. We normalized the gene expression by the total expression in each cell, multiplied it by 10 000 and log10-transformed it (hereafter, log-normalized expression level) using the NormalizeData function. We extracted genes with highly variable expression using the FindVariableGenes function and the following thresholds: mean of 0.1-3 and dispersion of log-normalized expression level of >0.5. We linear-regressed out the log-normalized expression level using the number of unique molecular identifiers in each cell, and scaled and centred residuals (hereafter, scaled expression level). We performed principal component analysis of the scaled expression level of the genes with highly variable expression and identified several PCs in each set of data that had a clearly larger standard deviation than the remaining PCs in PCElbowPlot function. We performed t-SNE using the RunTSNE function and the PCs with a large standard deviation as input and plotted the data for each cell in two dimensions. ChIP-seq analysis We reanalysed public data of ChIP-seq for Myrf (28). Briefly, in the originally reported work, cultured rat oligodendrocyte progenitor cells were transfected with vectors expressing either Myc-MYRF or untagged MYRF. After 48 h in differentiating conditions, ChIP-seq using Myc-tag antibody was performed, and 2102 and 17 peaks were detected by MACS (51) in Myc-MYRF and untagged MYRF experiments, respectively. Because all of the 17 peaks were detected among the 2102 peaks, there were 2085 peaks specific to the Myc-MYRF experiments. Using information on these 2085 peaks, we analysed the enrichment of biological process-related GO terms among the putative target genes of Myrf. For this purpose, we used CHIP-ENRICH (52) with the default settings. The Locus Definition parameter, which defines the threshold of distance between a Myrf peak and a TSS of a putative target gene, was determined as ‘<10 kb’, according to the relationship between the Myrf peaks and TSS. In this analysis, we analysed the enrichment of GO terms related to proliferation (number of terms, 14), migration (number of terms, 12) and EMT/MET (number of terms, 5), which were selected prior to the analysis and corrected the obtained raw P-values using the Benjamini–Hochberg method. To summarize the functions of all GO terms with a raw P-value of <0.05 in the CHIP-ENRICH results, we clustered the GO terms whose genes overlapped with >0.7 overlap coefficient using the EnrichmentMap (53) plugin of CytoScape v3.6.1 (54). To prevent the overclustering of unrelated GO terms, we excluded GO terms that contained >1000 genes. The distance from Myrf peaks to the nearest TSS outputted by the CHIP-ENRICH analysis was used to generate plots of the relationship between TSS and Myrf peaks. Acknowledgements We thank all of the participants for their cooperation in this research. We also thank Ms K. Takabe, Mr T. Miyama, Ms N. Watanabe, Ms M. Sato, Mr S. Nakamura and Ms S. Sugimoto at the Department of Human Genetics, Yokohama City University Graduate School of Medicine, for their technical assistance. We are also grateful to Edanz (www.edanzediting.co.jp) for editing the English text of a draft of this manuscript. Conflict of Interest statement. None declared. Funding Japan Agency for Medical Research and Development (AMED) (grant numbers JP18ek0109280, JP18dm0107090, JP18ek0109301, JP18ek0109348 and JP18kk020500 to N. Matsumoto); Japan Society for the Promotion of Science (JSPS) KAKENHI (grant numbers JP17H01539 to N. Matsumoto, JP16H05160 to H.S., JP16H05357 to N. Miyake, JP16H06254 to A.T., JP15K10367 to M.N., JP17K10080 to S.M., JP17K15630 to T.M., JP17H06994 to A.F.; Ministry of Health, Labour and Welfare (to N. Matsumoto); the Takeda Science Foundation (to N. Miyake, H.S, N. Matsumoto); and The Ichiro Kanehara Foundation for the Promotion of Medical Science & Medical Care (S.M.). References 1 Lee , P.A. , Houk , C.P. , Ahmed , S.F. and Hughes , I.A. ( 2006 ) Consensus statement on management of intersex disorders. International Consensus Conference on Intersex . Pediatrics , 118 , e488 – e500 . Google Scholar Crossref Search ADS PubMed WorldCat 2 Bashamboo , A. and McElreavey , K. ( 2013 ) Gene mutations associated with anomalies of human gonad formation . Sex. Dev. , 7 , 126 – 146 . Google Scholar Crossref Search ADS PubMed WorldCat 3 Dong , Y. , Yi , Y. , Yao , H. , Yang , Z. , Hu , H. , Liu , J. , Gao , C. , Zhang , M. , Zhou , L. , Asan , L. , Yi , X. and Liang , Z. ( 2016 ) Targeted next-generation sequencing identification of mutations in patients with disorders of sex development . BMC Med. Genet. , 17 , 23 . Google Scholar Crossref Search ADS PubMed WorldCat 4 Baxter , R.M. , Arboleda , V.A. , Lee , H. , Barseghyan , H. , Adam , M.P. , Fechner , P.Y. , Bargman , R. , Keegan , C. , Travers , S. , Schelley , S. et al. ( 2015 ) Exome sequencing for the diagnosis of 46,XY disorders of sex development . J. Clin. Endocrinol. Metab. , 100 , E333 – E344 . Google Scholar Crossref Search ADS PubMed WorldCat 5 Richards , S. , Aziz , N. , Bale , S. , Bick , D. , Das , S. , Gastier-Foster , J. , Grody , W.W. , Hegde , M. , Lyon , E. , Spector , E. , Voelkerding , K. and Rehm , H.L. ( 2015 ) Standards and guidelines for the interpretation of sequence variants: a joint consensus recommendation of the American College of Medical Genetics and Genomics and the Association for Molecular Pathology . Genet. Med. , 17 , 405 – 424 . Google Scholar Crossref Search ADS PubMed WorldCat 6 Samocha , K.E. , Robinson , E.B. , Sanders , S.J. , Stevens , C. , Sabo , A. , McGrath , L.M. , Kosmicki , J.A. , Rehnstrom , K. , Mallick , S. , Kirby , A. et al. ( 2014 ) A framework for the interpretation of de novo mutation in human disease . Nat. Genet. , 46 , 944 – 950 . Google Scholar Crossref Search ADS PubMed WorldCat 7 Adzhubei , I.A. ( 2010 ) A method and server for predicting damaging missense mutations . Nat. Methods , 7 , 248 – 249 . Google Scholar Crossref Search ADS PubMed WorldCat 8 Ng , P.C. and Henikoff , S. ( 2003 ) SIFT: predicting amino acid changes that affect protein function . Nucleic Acids Res. , 31 , 3812 – 3814 . Google Scholar Crossref Search ADS PubMed WorldCat 9 Schwarz , J.M. , Cooper , D.N. , Schuelke , M. and Seelow , D. ( 2014 ) MutationTaster2: mutation prediction for the deep-sequencing age . Nat. Methods , 11 , 361 – 362 . Google Scholar Crossref Search ADS PubMed WorldCat 10 Kircher , M. , Witten , D.M. , Jain , P. , O'Roak , B.J. , Cooper , G.M. and Shendure , J. ( 2014 ) A general framework for estimating the relative pathogenicity of human genetic variants . Nat. Genet. , 46 , 310 – 315 . Google Scholar Crossref Search ADS PubMed WorldCat 11 Zhen , X. , Li , B. , Hu , F. , Yan , S. , Meloni , G. , Li , H. and Shi , N. ( 2017 ) Crystal structure of the DNA-binding domain of myelin-gene regulatory factor . Sci. Rep. , 7 , 3696 . Google Scholar Crossref Search ADS PubMed WorldCat 12 Firth , H.V. , Richards , S.M. , Bevan , A.P. , Clayton , S. , Corpas , M. , Rajan , D. , Van Vooren , S. , Moreau , Y. , Pettett , R.M. and Carter , N.P. ( 2009 ) DECIPHER: database of chromosomal imbalance and phenotype in humans using Ensembl resources . Am. J. Hum. Genet. , 84 , 524 – 533 . Google Scholar Crossref Search ADS PubMed WorldCat 13 Kardon , G. , Ackerman , K.G. , McCulley , D.J. , Shen , Y. , Wynn , J. , Shang , L. , Bogenschutz , E. , Sun , X. and Chung , W.K. ( 2017 ) Congenital diaphragmatic hernias: from genes to mechanisms to therapies . Dis. Model. Mech. , 10 , 955 – 970 . Google Scholar Crossref Search ADS PubMed WorldCat 14 Orvis , G.D. and Behringer , R.R. ( 2007 ) Cellular mechanisms of Mullerian duct formation in the mouse . Dev. Biol. , 306 , 493 – 504 . Google Scholar Crossref Search ADS PubMed WorldCat 15 Guioli , S. , Sekido , R. and Lovell-Badge , R. ( 2007 ) The origin of the Mullerian duct in chick and mouse . Dev. Biol. , 302 , 389 – 398 . Google Scholar Crossref Search ADS PubMed WorldCat 16 Carmona , R. , Canete , A. , Cano , E. , Ariza , L. , Rojas , A. and Munoz-Chapuli , R. ( 2016 ) Conditional deletion of WT1 in the septum transversum mesenchyme causes congenital diaphragmatic hernia in mice . Elife , 5 . doi: https://doi.org/10.7554/eLife.16009 . WorldCat 17 Piprek , R.P. ( 2016 ) Molecular Mechanisms of Cell Differentiation in Gonad Development . Springer Publishing , NY . 18 Ariza , L. , Carmona , R. , Canete , A. , Cano , E. and Munoz-Chapuli , R. ( 2016 ) Coelomic epithelium-derived cells in visceral morphogenesis . Dev. Dyn. , 245 , 307 – 322 . Google Scholar Crossref Search ADS PubMed WorldCat 19 Li , L. , Dong , J. , Yan , L. , Yong , J. , Liu , X. , Hu , Y. , Fan , X. , Wu , X. , Guo , H. , Wang , X. et al. ( 2017 ) Single-cell RNA-seq analysis maps development of human germline cells and gonadal niche interactions . Cell Stem Cell , 20 , 858 – 873 e854 . Google Scholar Crossref Search ADS PubMed WorldCat 20 Rudat , C. , Grieskamp , T. , Rohr , C. , Airik , R. , Wrede , C. , Hegermann , J. , Herrmann , B.G. , Schuster-Gossler , K. and Kispert , A. ( 2014 ) Upk3b is dispensable for development and integrity of urothelium and mesothelium . PLoS One , 9 , e112112. doi: https://doi.org/10.1371/journal.pone.0112112 . WorldCat 21 Mork , L. , Maatouk , D.M. , McMahon , J.A. , Guo , J.J. , Zhang , P. , McMahon , A.P. and Capel , B. ( 2012 ) Temporal differences in granulosa cell specification in the ovary reflect distinct follicle fates in mice . Biol. Reprod. , 86 , 37 . Google Scholar Crossref Search ADS PubMed WorldCat 22 Rastetter , R.H. , Bernard , P. , Palmer , J.S. , Chassot , A.A. , Chen , H. , Western , P.S. , Ramsay , R.G. , Chaboissier , M.C. and Wilhelm , D. ( 2014 ) Marker genes identify three somatic cell types in the fetal mouse ovary . Dev. Biol. , 394 , 242 – 252 . Google Scholar Crossref Search ADS PubMed WorldCat 23 Lee , H.J. , Pazin , D.E. , Kahlon , R.S. , Correa , S.M. and Albrecht , K.H. ( 2009 ) Novel markers of early ovarian pre-granulosa cells are expressed in an Sry-like pattern . Dev. Dyn. , 238 , 812 – 825 . Google Scholar Crossref Search ADS PubMed WorldCat 24 Stockinger , H. , Gadd , S.J. , Eher , R. , Majdic , O. , Schreiber , W. , Kasinrerk , W. , Strass , B. , Schnabl , E. and Knapp , W. ( 1990 ) Molecular characterization and functional analysis of the leukocyte surface protein CD31 . J. Immunol. , 145 , 3889 – 3897 . Google Scholar PubMed WorldCat 25 Gory , S. , Vernet , M. , Laurent , M. , Dejana , E. , Dalmon , J. and Huber , P. ( 1999 ) The vascular endothelial-cadherin promoter directs endothelial-specific expression in transgenic mice . Blood , 93 , 184 – 192 . Google Scholar PubMed WorldCat 26 de Santa Barbara , P. , Moniot , B. , Poulat , F. and Berta , P. ( 2000 ) Expression and subcellular localization of SF-1, SOX9, WT1, and AMH proteins during early human testicular development . Dev. Dyn. , 217 , 293 – 298 . Google Scholar Crossref Search ADS PubMed WorldCat 27 van den Driesche , S. , Walker , M. , McKinnell , C. , Scott , H.M. , Eddie , S.L. , Mitchell , R.T. , Seckl , J.R. , Drake , A.J. , Smith , L.B. , Anderson , R.A. and Sharpe , R.M. ( 2012 ) Proposed role for COUP-TFII in regulating fetal Leydig cell steroidogenesis, perturbation of which leads to masculinization disorders in rodents . PLoS One , 7 , e37064. doi: https://doi.org/10.1371/journal.pone.0037064 . WorldCat 28 Miyabayashi , K. , Katoh-Fukui , Y. , Ogawa , H. , Baba , T. , Shima , Y. , Sugiyama , N. , Kitamura , K. and Morohashi , K. ( 2013 ) Aristaless related homeobox gene, Arx, is implicated in mouse fetal Leydig cell differentiation possibly through expressing in the progenitor cells . PLoS One , 8 , e68050. doi: https://doi.org/10.1371/journal.pone.0068050 . WorldCat 29 Miyabayashi , K. , Shima , Y. , Inoue , M. , Sato , T. , Baba , T. , Ohkawa , Y. , Suyama , M. and Morohashi , K.I. ( 2017 ) Alterations in fetal Leydig cell gene expression during fetal and adult development . Sex. Dev. , 11 , 53 – 63 . Google Scholar Crossref Search ADS PubMed WorldCat 30 Han , X. , Wang , R. , Zhou , Y. , Fei , L. , Sun , H. , Lai , S. , Saadatpour , A. , Zhou , Z. , Chen , H. , Ye , F. et al. ( 2018 ) Mapping the mouse cell atlas by microwell-seq . Cell , 172 , 1091 – 1107 . Google Scholar Crossref Search ADS PubMed WorldCat 31 Bujalka , H. , Koenning , M. , Jackson , S. , Perreau , V.M. , Pope , B. , Hay , C.M. , Mitew , S. , Hill , A.F. , Lu , Q.R. , Wegner , M. et al. ( 2013 ) MYRF is a membrane-associated transcription factor that autoproteolytically cleaves to directly activate myelin genes . PLoS Biol , 11 , e1001625. doi: https://doi.org/10.1371/journal.pbio.1001625 . WorldCat 32 Emery , B. , Agalliu , D. , Cahoy , J.D. , Watkins , T.A. , Dugas , J.C. , Mulinyawe , S.B. , Ibrahim , A. , Ligon , K.L. , Rowitch , D.H. and Barres , B.A. ( 2009 ) Myelin gene regulatory factor is a critical transcriptional regulator required for CNS myelination . Cell , 138 , 172 – 185 . Google Scholar Crossref Search ADS PubMed WorldCat 33 Handstad , T. , Rye , M. , Mocnik , R. , Drablos , F. and Saetrom , P. ( 2012 ) Cell-type specificity of ChIP-predicted transcription factor binding sites . BMC Genomics , 13 , 372 . Google Scholar Crossref Search ADS PubMed WorldCat 34 Pinz , H. , Pyle , L.C. , Li , D. , Izumi , K. , Skraban , C. , Tarpinian , J. , Braddock , S.R. , Telegrafi , A. , Monaghan , K.G. , Zackai , E. and Bhoj , E.J. ( 2018 ) De novo variants in myelin regulatory factor (MYRF) as candidates of a new syndrome of cardiac and urogenital anomalies . Am. J. Med. Genet. A , 176 , 969 – 972 . Google Scholar Crossref Search ADS PubMed WorldCat 35 Sarto , G.E. and Opitz , J.M. ( 1973 ) The XY gonadal agenesis syndrome . J. Med. Genet. , 10 , 288 – 293 . Google Scholar Crossref Search ADS PubMed WorldCat 36 Levinson , G. , Zarate , A. , Guzman-Toledano , R. , Canales , E.S. and Jimenez , M. ( 1976 ) An XX female with sexual infantilism, absent gonads, and lack of Mullerian ducts . J. Med. Genet. , 13 , 68 – 69 . Google Scholar Crossref Search ADS PubMed WorldCat 37 De Leon , F.D. , Hersh , J.H. , Sanfilippo , J.S. , Schikler , K.N. and Yen , F.F. ( 1984 ) Gonadal and Mullerian duct agenesis in a girl with 46,X,i(Xq) . Obstet. Gynecol. , 63 , 81s – 83s . Google Scholar PubMed WorldCat 38 Koenning , M. , Jackson , S. , Hay , C.M. , Faux , C. , Kilpatrick , T.J. , Willingham , M. and Emery , B. ( 2012 ) Myelin gene regulatory factor is required for maintenance of myelin and mature oligodendrocyte identity in the adult CNS . J. Neurosci. , 32 , 12528 – 12542 . Google Scholar Crossref Search ADS PubMed WorldCat 39 Kurahashi , H. , Azuma , Y. , Masuda , A. , Okuno , T. , Nakahara , E. , Imamura , T. , Saitoh , M. , Mizuguchi , M. , Shimizu , T. , Ohno , K. and Okumura , A. ( 2018 ) MYRF is associated with encephalopathy with reversible myelin vacuolization . Ann. Neurol. , 83 , 98 – 106 . Google Scholar Crossref Search ADS PubMed WorldCat 40 Hamanaka , K. , Miyatake , S. , Zerem , A. , Lev , D. , Blumkin , L. , Yokochi , K. , Fujita , A. , Imagawa , E. , Iwama , K. , Nakashima , M. et al. ( 2018 ) Expanding the phenotype of IBA57 mutations: related leukodystrophy can remain asymptomatic . J. Hum. Genet. , 63 , 1223 – 1229 . Google Scholar Crossref Search ADS PubMed WorldCat 41 McKenna , A. , Hanna , M. , Banks , E. , Sivachenko , A. , Cibulskis , K. , Kernytsky , A. , Garimella , K. , Altshuler , D. , Gabriel , S. , Daly , M. and DePristo , M.A. ( 2010 ) The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data . Genome Res. , 20 , 1297 – 1303 . Google Scholar Crossref Search ADS PubMed WorldCat 42 DePristo , M.A. , Banks , E. , Poplin , R. , Garimella , K.V. , Maguire , J.R. , Hartl , C. , Philippakis , A.A. , del Angel , G. , Rivas , M.A. , Hanna , M. et al. ( 2011 ) A framework for variation discovery and genotyping using next-generation DNA sequencing data . Nat. Genet. , 43 , 491 – 498 . Google Scholar Crossref Search ADS PubMed WorldCat 43 Wang , K. , Li , M. and Hakonarson , H. ( 2010 ) ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data . Nucleic Acids Res. , 38 , e164 . Google Scholar Crossref Search ADS PubMed WorldCat 44 Cingolani , P. , Platts , A. , Wang le , L. , Coon , M. , Nguyen , T. , Wang , L. , Land , S.J. , Lu , X. and Ruden , D.M. ( 2012 ) A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3 . Fly (Austin) , 6 , 80 – 92 . Google Scholar Crossref Search ADS PubMed WorldCat 45 Liu , X. , Jian , X. and Boerwinkle , E. ( 2011 ) dbNSFP: a lightweight database of human nonsynonymous SNPs and their functional predictions . Hum. Mutat. , 32 , 894 – 899 . Google Scholar Crossref Search ADS PubMed WorldCat 46 Lek , M. ( 2016 ) Analysis of protein-coding genetic variation in 60,706 humans . Nature , 536 , 285 – 291 . Google Scholar Crossref Search ADS PubMed WorldCat 47 Nagasaki , M. , Yasuda , J. , Katsuoka , F. , Nariai , N. , Kojima , K. , Kawai , Y. , Yamaguchi-Kabata , Y. , Yokozawa , J. , Danjoh , I. , Saito , S. et al. ( 2015 ) Rare variant discovery by deep whole-genome sequencing of 1,070 Japanese individuals . Nat. Commun. , 6 , 8018 . Google Scholar Crossref Search ADS PubMed WorldCat 48 Li , B. and Leal , S.M. ( 2008 ) Methods for detecting associations with rare variants for common diseases: application to analysis of sequence data . Am. J. Hum. Genet. , 83 , 311 – 321 . Google Scholar Crossref Search ADS PubMed WorldCat 49 Zhan , X. , Hu , Y. , Li , B. , Abecasis , G.R. and Liu , D.J. ( 2016 ) RVTESTS: an efficient and comprehensive tool for rare variant association analysis using sequence data . Bioinformatics , 32 , 1423 – 1426 . Google Scholar Crossref Search ADS PubMed WorldCat 50 Butler , A. , Hoffman , P. , Smibert , P. , Papalexi , E. and Satija , R. ( 2018 ) Integrating single-cell transcriptomic data across different conditions, technologies, and species . Nat. Biotechnol. , 36 , 411 – 420 . Google Scholar Crossref Search ADS PubMed WorldCat 51 Zhang , Y. , Liu , T. , Meyer , C.A. , Eeckhoute , J. , Johnson , D.S. , Bernstein , B.E. , Nusbaum , C. , Myers , R.M. , Brown , M. , Li , W. and Liu , X.S. ( 2008 ) Model-based analysis of ChIP-seq (MACS) . Genome Biol. , 9 , R137 . Google Scholar Crossref Search ADS PubMed WorldCat 52 Welch , R.P. , Lee , C. , Imbriano , P.M. , Patil , S. , Weymouth , T.E. , Smith , R.A. , Scott , L.J. and Sartor , M.A. ( 2014 ) ChIP-Enrich: gene set enrichment testing for ChIP-seq data . Nucleic Acids Res. , 42 , e105 . Google Scholar Crossref Search ADS PubMed WorldCat 53 Merico , D. , Isserlin , R. , Stueker , O. , Emili , A. and Bader , G.D. ( 2010 ) Enrichment map: a network-based method for gene-set enrichment visualization and interpretation . PLoS One , 5 , e13984. WorldCat 54 Shannon , P. , Markiel , A. , Ozier , O. , Baliga , N.S. , Wang , J.T. , Ramage , D. , Amin , N. , Schwikowski , B. and Ideker , T. ( 2003 ) Cytoscape: a software environment for integrated models of biomolecular interaction networks . Genome Res. , 13 , 2498 – 2504 . Google Scholar Crossref Search ADS PubMed WorldCat © The Author(s) 2019. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model) http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Human Molecular Genetics Oxford University Press

MYRF haploinsufficiency causes 46,XY and 46,XX disorders of sex development: bioinformatics consideration

Download PDF
11 pages

Loading next page...
 
/lp/oxford-university-press/myrf-haploinsufficiency-causes-46-xy-and-46-xx-disorders-of-sex-GQNPvmVacL

References (56)

Publisher
Oxford University Press
Copyright
© The Author(s) 2019. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oup.com
ISSN
0964-6906
eISSN
1460-2083
DOI
10.1093/hmg/ddz066
Publisher site
See Article on Publisher Site

Abstract

Abstract Disorders of sex development (DSDs) are defined as congenital conditions in which chromosomal, gonadal or anatomical sex is atypical. In many DSD cases, genetic causes remain to be elucidated. Here, we performed a case–control exome sequencing study comparing gene-based burdens of rare damaging variants between 26 DSD cases and 2625 controls. We found exome-wide significant enrichment of rare heterozygous truncating variants in the MYRF gene encoding myelin regulatory factor, a transcription factor essential for oligodendrocyte development. All three variants occurred de novo. We identified an additional 46,XY DSD case of a de novo damaging missense variant in an independent cohort. The clinical symptoms included hypoplasia of Müllerian derivatives and ovaries in 46,XX DSD patients, defective development of Sertoli and Leydig cells in 46,XY DSD patients and congenital diaphragmatic hernia in one 46,XY DSD patient. As all of these cells and tissues are or partly consist of coelomic epithelium (CE)-derived cells (CEDC) and CEDC developed from CE via proliferaiton and migration, MYRF might be related to these processes. Consistent with this hypothesis, single-cell RNA sequencing of foetal gonads revealed high expression of MYRF in CE and CEDC. Reanalysis of public chromatin immunoprecipitation sequencing data for rat Myrf showed that genes regulating proliferation and migration were enriched among putative target genes of Myrf. These results suggested that MYRF is a novel causative gene of 46,XY and 46,XX DSD and MYRF is a transcription factor regulating CD and/or CEDC proliferation and migration, which is essential for development of multiple organs. Introduction Disorders of sex development (DSDs) are defined as congenital conditions in which chromosomal, gonadal or anatomical sex is atypical (1). As per this definition, sex chromosome abnormalities are a type of DSD and include Klinefelter syndrome (47,XXY) with infertility, Turner syndrome (45,XO) with streak gonads, 45,X/46,XY DSD with nearly normal male to Turner syndrome-like phenotype and 46,XX/46,XY DSD with ovotestis. Other types of DSD, namely, 46,XY and 46,XX DSD, are known to be caused by pathogenic variants in specific genes including NR5A1, SRY and SOX9 (2). Recent targeted next-generation sequencing of known mutated genes or whole-exome sequencing (WES) has reported diagnostic yields of 35–40% in DSD cases without sex chromosome aneuploidy (3,4). Thus, in the large proportion of sex chromosome aneuploidy-negative DSD cases, the genetic cause remains unclear. Figure 1 View largeDownload slide Enrichment of rare truncating variants of MYRF in DSD cases. (A, B) Statistical significance (P-value) of gene-based burden of rare damaging variants in DSD cases compared with controls. Two types of damaging variants, LGD model (A) and LGD + D-mis model (B), are shown. Manhattan plot (left) and QQ plot (right) of P-value are shown. Red lines in Manhattan plots, thresholds for statistical significance; red lines in QQ plots, x = y. (C) Four pedigrees with DSD cases harbouring MYRF damaging variants.The properties of the MYRF variants are shown below each individual. Arrows, probands; squares, males; circles, females; closed symbols, affected; open symbols, unaffected; wt, wild-type allele. (D) Poisson distribution of the number of de novo truncating variants of MYRF in 26 samples under the null model of mutation (6). The probability of observing three or more de novo truncating variants of MYRF in 26 individuals (52 alleles) was 1.95E-012. (E) Evolutionary conservation of MYRF Gln443 that is mutated in Patient 5. Note that Gln443 was near an interface for trimerization (11). (F) Hydrogen bonds formed by p.Gln443Pro. The side chain of wild-type Gln443 formed hydrogen bonds with Lys437 and Ala440 (upper), while that of p.Gln443Pro formed two novel hydrogen bonds (pink dotted line) with Ala440. Figure 1 View largeDownload slide Enrichment of rare truncating variants of MYRF in DSD cases. (A, B) Statistical significance (P-value) of gene-based burden of rare damaging variants in DSD cases compared with controls. Two types of damaging variants, LGD model (A) and LGD + D-mis model (B), are shown. Manhattan plot (left) and QQ plot (right) of P-value are shown. Red lines in Manhattan plots, thresholds for statistical significance; red lines in QQ plots, x = y. (C) Four pedigrees with DSD cases harbouring MYRF damaging variants.The properties of the MYRF variants are shown below each individual. Arrows, probands; squares, males; circles, females; closed symbols, affected; open symbols, unaffected; wt, wild-type allele. (D) Poisson distribution of the number of de novo truncating variants of MYRF in 26 samples under the null model of mutation (6). The probability of observing three or more de novo truncating variants of MYRF in 26 individuals (52 alleles) was 1.95E-012. (E) Evolutionary conservation of MYRF Gln443 that is mutated in Patient 5. Note that Gln443 was near an interface for trimerization (11). (F) Hydrogen bonds formed by p.Gln443Pro. The side chain of wild-type Gln443 formed hydrogen bonds with Lys437 and Ala440 (upper), while that of p.Gln443Pro formed two novel hydrogen bonds (pink dotted line) with Ala440. To identify a novel causative gene of DSD, we here performed a case–control WES study comparing the gene-based burden of rare damaging variants between DSD cases (N = 26) and healthy controls (N = 2625). We identified MYRF (MIM: 608329), encoding myelin regulatory factor, as a causative gene of DSD with exome-wide statistical significance. We consequently show convergent evidence indicating a critical role of MYRF in the development of coelomic epithelium (CE)-derived cells (CEDC), which is supported by detailed phenotypic assessment of the DSD cases with the MYRF variants, single-cell gene expression analysis and comprehensive mapping of putative target genes of the transcription factor MYRF. Result In this study, we performed WES on 26 DSD cases and 2625 controls. For the DSD cases, we confirmed that none of the individuals harboured known genetic causes, namely, sex chromosome aneuploidy (Supplementary Material, Fig. S1A and B) or pathogenic/likely pathogenic variants in 68 human-curated genes implicated in DSD (Supplementary Material, Table S1) (2,4). The classification of variant pathogenicity was based on the guidelines from the American College of Medical Genetics and Genomics (5). By using the obtained WES data, we confirmed that there was no apparent population stratification (Supplementary Material, Fig. S1C), no strong relatedness among samples (Supplementary Material, Fig. S1D) and no sample contamination (Supplementary Material, Fig. S1E) in any of the DSD cases and controls, confirming that the quality of our data was sufficient for the following case–control study. We next performed an exome-wide analysis comparing the gene-based burden of rare damaging variants between cases and controls (see Materials and Methods). We analysed two types of damaging variants: likely gene-disrupting variants (LGDs) only or LGD + damaging missense variants (D-mis). For LGD burden, one gene, MYRF at 11q12.2 showed a P-value (8.38E−07) below the threshold for exome-wide statistical significance (2.05E-06, see Materials and Methods) (Fig. 1A, left). For MYRF, there were three LGD variants in cases and none in controls. A quantile–quantile (QQ) plot of observed and expected P-values revealed that there was no inflation in the observed P-values and identified MYRF as the only gene with upward deviation from the theoretical uniform distribution of P-values (Fig. 1A, right). For LGD + D-mis burden, no gene had a P-value below the threshold described above (Fig. 1B, left), and there was no upward outlier from the theoretical distribution of P-value (Fig. 1B, right). Three LGD variants were identified in MYRF (NM_001127392.2): c.278del, p.(Pro93Argfs*7); c.789del, p.(Ser264Alafs*8); and c.1642_1666del, p.(Ala548Thrfs*49) at the fourth, fifth and eighth exons, respectively. All of these variants not located at the last exon or the last 50 bp from the penultimate exon are expected to be subjected to nonsense-mediated mRNA decay, although we could not experimentally test the effect of these LGD variants because cell samples were unavailable in patients. We confirmed the existence of these variants by Sanger sequencing and revealed that all of them are de novo variants (Fig. 1C). In Family 2, the c.789del: p.(Ser264Alafs*8) variant was shared with the affected monozygotic twin. In Family 3, the c.1642_1666del: p.(Ala548Thrfs*49) variant was not found in two unaffected siblings (Fig. 1C). Because all of the LGD variants in MYRF occurred de novo, we additionally investigated whether the observed number of de novo MYRF variants was significantly more than expected in 26 DSD cases under the null model of mutation (6). Under this model, the probability of observing three or more de novo LGD variants in MYRF in 26 individuals was 1.95E-12 (Fig. 1D) and below the threshold of statistical significance (2.05E-06, see Materials and Methods). This result further demonstrates the robust statistical association between LGD variants in MYRF and DSD. The statistical evidence for the enrichment of de novo MYRF variants in DSD motivated us to screen rare MYRF variants in another DSD cohort. Of 68 DSD cases of unknown genetic cause, we identified a de novo rare MYRF missense variant, c.1328A>C: p.(Gln443Pro), in Family 4 (Fig. 1C). The pathogenicity of this variant was supported by multiple lines of evidence as follows. (1) Most in silico tools predicted the variant to be damaging (Supplementary Material, Table S2) (7–10). (2) The amino acid at the variant was conserved among species (Fig. 1E). (3) There is no missense variant affecting the same amino acid (Gln443) in general populations: ESP6500, gnomAD and ToMMo3.5KJPN. (4) The variant was located adjacent to the interface for homo-trimerization (11) (Fig. 1E), where missense variants are depleted in general populations, as shown in the DECIPHER browser (12) (data not shown). The missense variant was predicted by Swiss PDB viewer structural analysis (Protein Data Bank ID: 5H5P) to intramolecularly form a novel hydrogen bond (Fig. 1F) and might change the structure of the adjacent interface for homo-trimerization. When we analysed the probability of observing four or more de novo LGD or missense variants in MYRF among 94 (26 + 68) individuals under the null model of mutation (6), we calculated the probability as 1.0E-10, which was again below the threshold of exome-wide statistical significance (2.05E-06). Table 1 Summary of clinical and genetic findings in DSD cases with MYRF damaging variants Case Pt1 Monozygotic twin Pt4 Pt5 Pt2 Pt3  Karyotype 46,XY 46,XX 46,XX 46,XY 46,XY  Social sex F—>M at 4 y F F M M  Present age 14 years 21 years 21 years 4 years 4 months 2 years 5 months  Genital findings   Age at examination 2 months 15 years 15 years 3 momths 0 (At birth)   Tanner stage I III (breast); II (pubis) II (breast); II (pubis) Prepubertal Prepubertal   Testis size (mL) Unknown - - Small (<1 mL, L) 0.9 × 0.5 × 0.6 cm on echo Small (<1 mL, L and R)   Cryptorchidism Abdominal (R); inguinal (L) - - Abdominal (L and R) Inguinal (L and R)   Vanishing testis No - - + (R) No   Hypospadias Perineal type - - Penile type Penoscrotal type   Microphallus 2.3 cm (−3.0 SD) - - 1.8 cm (−3.0 SD) 1.5 cm (<−3.0 SD)   Abnormal scrotum Anterior-positioned and bifid - - No Small and bifid   Ovary - Absent (L); small (R) Absent (L); small (R) - -   Fallopian tube - Absent Absent - -   Uterus - Absent (L); restiform cervix (R) Absent (L); restiform cervix (R) - -   Vagina - Absent, surrounded by a cyst Absent, surrounded by a cyst - -  Endocrine findings   Age at examination 2 months 15 years 15 years 1 year 0 month 3 months 7 months   LH (IU/L) baseline NA 45.2 25.7 0.73 (0.2–1.9) 8.9 (0.2–1.9) 0.3 (0.2–1.9)   LH (IU/L) GnRH-stimulated NA NA NA 9.49 (1.1–6.0) NA 8.4 (1.1–6.0)   FSH (IU/L) baseline NA 70.1 113.1 3.91 (<0.3–2.4) 8.9 (<0.3–2.4) 3.4 (<0.3–2.4)   FSH (IU/L) GnRH-stimulated NA NA NA 21.62 (1.9–7.6) NA 23.9 (1.9–7.6)   T (ng/dL) baseline 125 (120 ~ 400) 55.0 54.9 <5 (<5) 116 (120–400) <5 (<5)   T (ng/dL) hCG-stimulated 440 (> 200) NA NA 518 (> 200) NA 78 (>200)   DHT (ng/dL) baseline NA NA NA <5 (<5) 20 NA   AMH (ng/mL) baseline NA NA NA 18.3 (43.3–79.3) 19.6 (43.3–79.3) 15.0 (43.3–79.3)  Treatment   Orchidopexy (age) 1 years - - 1 year 7 months Not yet   Genitoplasty (age) 2 years 15 years 15 years 2 years 11 month Not yet   TE (25 mg i.m.) (age) 3 × (6 ~ 10 months) - - No 1 × (5 months)   DHT cream (age) No - - Yes (2 years 8 months) No   PL increment (cm) 2.3→3.9 (+ 1.1 SD) - - Increased (no record) 1.5→3.0 (+1.2 SD)  Central nervous system   Developmental delay No No No Yes No   DQ NA NA NA 63 NA   Seizure No No No No No   ASD/ADHD No No No NA No  Others   Diaphragmatic hernia No No No Yes (L) No Case Pt1 Monozygotic twin Pt4 Pt5 Pt2 Pt3  Karyotype 46,XY 46,XX 46,XX 46,XY 46,XY  Social sex F—>M at 4 y F F M M  Present age 14 years 21 years 21 years 4 years 4 months 2 years 5 months  Genital findings   Age at examination 2 months 15 years 15 years 3 momths 0 (At birth)   Tanner stage I III (breast); II (pubis) II (breast); II (pubis) Prepubertal Prepubertal   Testis size (mL) Unknown - - Small (<1 mL, L) 0.9 × 0.5 × 0.6 cm on echo Small (<1 mL, L and R)   Cryptorchidism Abdominal (R); inguinal (L) - - Abdominal (L and R) Inguinal (L and R)   Vanishing testis No - - + (R) No   Hypospadias Perineal type - - Penile type Penoscrotal type   Microphallus 2.3 cm (−3.0 SD) - - 1.8 cm (−3.0 SD) 1.5 cm (<−3.0 SD)   Abnormal scrotum Anterior-positioned and bifid - - No Small and bifid   Ovary - Absent (L); small (R) Absent (L); small (R) - -   Fallopian tube - Absent Absent - -   Uterus - Absent (L); restiform cervix (R) Absent (L); restiform cervix (R) - -   Vagina - Absent, surrounded by a cyst Absent, surrounded by a cyst - -  Endocrine findings   Age at examination 2 months 15 years 15 years 1 year 0 month 3 months 7 months   LH (IU/L) baseline NA 45.2 25.7 0.73 (0.2–1.9) 8.9 (0.2–1.9) 0.3 (0.2–1.9)   LH (IU/L) GnRH-stimulated NA NA NA 9.49 (1.1–6.0) NA 8.4 (1.1–6.0)   FSH (IU/L) baseline NA 70.1 113.1 3.91 (<0.3–2.4) 8.9 (<0.3–2.4) 3.4 (<0.3–2.4)   FSH (IU/L) GnRH-stimulated NA NA NA 21.62 (1.9–7.6) NA 23.9 (1.9–7.6)   T (ng/dL) baseline 125 (120 ~ 400) 55.0 54.9 <5 (<5) 116 (120–400) <5 (<5)   T (ng/dL) hCG-stimulated 440 (> 200) NA NA 518 (> 200) NA 78 (>200)   DHT (ng/dL) baseline NA NA NA <5 (<5) 20 NA   AMH (ng/mL) baseline NA NA NA 18.3 (43.3–79.3) 19.6 (43.3–79.3) 15.0 (43.3–79.3)  Treatment   Orchidopexy (age) 1 years - - 1 year 7 months Not yet   Genitoplasty (age) 2 years 15 years 15 years 2 years 11 month Not yet   TE (25 mg i.m.) (age) 3 × (6 ~ 10 months) - - No 1 × (5 months)   DHT cream (age) No - - Yes (2 years 8 months) No   PL increment (cm) 2.3→3.9 (+ 1.1 SD) - - Increased (no record) 1.5→3.0 (+1.2 SD)  Central nervous system   Developmental delay No No No Yes No   DQ NA NA NA 63 NA   Seizure No No No No No   ASD/ADHD No No No NA No  Others   Diaphragmatic hernia No No No Yes (L) No NOTE: LH indicates luteinizing hormone; i.m., intramuscular injection; FSH, follicle-stimulating hormone; DHT, dihydrotestosterone; TE, testosterone enanthate; PL, penile length; DQ, development quotient measured with Kyoto Scale of Psychological Development; ASD, autism spectrum disorder; ADHD, attention-deficit hyperactivity disorder. View Large Table 1 Summary of clinical and genetic findings in DSD cases with MYRF damaging variants Case Pt1 Monozygotic twin Pt4 Pt5 Pt2 Pt3  Karyotype 46,XY 46,XX 46,XX 46,XY 46,XY  Social sex F—>M at 4 y F F M M  Present age 14 years 21 years 21 years 4 years 4 months 2 years 5 months  Genital findings   Age at examination 2 months 15 years 15 years 3 momths 0 (At birth)   Tanner stage I III (breast); II (pubis) II (breast); II (pubis) Prepubertal Prepubertal   Testis size (mL) Unknown - - Small (<1 mL, L) 0.9 × 0.5 × 0.6 cm on echo Small (<1 mL, L and R)   Cryptorchidism Abdominal (R); inguinal (L) - - Abdominal (L and R) Inguinal (L and R)   Vanishing testis No - - + (R) No   Hypospadias Perineal type - - Penile type Penoscrotal type   Microphallus 2.3 cm (−3.0 SD) - - 1.8 cm (−3.0 SD) 1.5 cm (<−3.0 SD)   Abnormal scrotum Anterior-positioned and bifid - - No Small and bifid   Ovary - Absent (L); small (R) Absent (L); small (R) - -   Fallopian tube - Absent Absent - -   Uterus - Absent (L); restiform cervix (R) Absent (L); restiform cervix (R) - -   Vagina - Absent, surrounded by a cyst Absent, surrounded by a cyst - -  Endocrine findings   Age at examination 2 months 15 years 15 years 1 year 0 month 3 months 7 months   LH (IU/L) baseline NA 45.2 25.7 0.73 (0.2–1.9) 8.9 (0.2–1.9) 0.3 (0.2–1.9)   LH (IU/L) GnRH-stimulated NA NA NA 9.49 (1.1–6.0) NA 8.4 (1.1–6.0)   FSH (IU/L) baseline NA 70.1 113.1 3.91 (<0.3–2.4) 8.9 (<0.3–2.4) 3.4 (<0.3–2.4)   FSH (IU/L) GnRH-stimulated NA NA NA 21.62 (1.9–7.6) NA 23.9 (1.9–7.6)   T (ng/dL) baseline 125 (120 ~ 400) 55.0 54.9 <5 (<5) 116 (120–400) <5 (<5)   T (ng/dL) hCG-stimulated 440 (> 200) NA NA 518 (> 200) NA 78 (>200)   DHT (ng/dL) baseline NA NA NA <5 (<5) 20 NA   AMH (ng/mL) baseline NA NA NA 18.3 (43.3–79.3) 19.6 (43.3–79.3) 15.0 (43.3–79.3)  Treatment   Orchidopexy (age) 1 years - - 1 year 7 months Not yet   Genitoplasty (age) 2 years 15 years 15 years 2 years 11 month Not yet   TE (25 mg i.m.) (age) 3 × (6 ~ 10 months) - - No 1 × (5 months)   DHT cream (age) No - - Yes (2 years 8 months) No   PL increment (cm) 2.3→3.9 (+ 1.1 SD) - - Increased (no record) 1.5→3.0 (+1.2 SD)  Central nervous system   Developmental delay No No No Yes No   DQ NA NA NA 63 NA   Seizure No No No No No   ASD/ADHD No No No NA No  Others   Diaphragmatic hernia No No No Yes (L) No Case Pt1 Monozygotic twin Pt4 Pt5 Pt2 Pt3  Karyotype 46,XY 46,XX 46,XX 46,XY 46,XY  Social sex F—>M at 4 y F F M M  Present age 14 years 21 years 21 years 4 years 4 months 2 years 5 months  Genital findings   Age at examination 2 months 15 years 15 years 3 momths 0 (At birth)   Tanner stage I III (breast); II (pubis) II (breast); II (pubis) Prepubertal Prepubertal   Testis size (mL) Unknown - - Small (<1 mL, L) 0.9 × 0.5 × 0.6 cm on echo Small (<1 mL, L and R)   Cryptorchidism Abdominal (R); inguinal (L) - - Abdominal (L and R) Inguinal (L and R)   Vanishing testis No - - + (R) No   Hypospadias Perineal type - - Penile type Penoscrotal type   Microphallus 2.3 cm (−3.0 SD) - - 1.8 cm (−3.0 SD) 1.5 cm (<−3.0 SD)   Abnormal scrotum Anterior-positioned and bifid - - No Small and bifid   Ovary - Absent (L); small (R) Absent (L); small (R) - -   Fallopian tube - Absent Absent - -   Uterus - Absent (L); restiform cervix (R) Absent (L); restiform cervix (R) - -   Vagina - Absent, surrounded by a cyst Absent, surrounded by a cyst - -  Endocrine findings   Age at examination 2 months 15 years 15 years 1 year 0 month 3 months 7 months   LH (IU/L) baseline NA 45.2 25.7 0.73 (0.2–1.9) 8.9 (0.2–1.9) 0.3 (0.2–1.9)   LH (IU/L) GnRH-stimulated NA NA NA 9.49 (1.1–6.0) NA 8.4 (1.1–6.0)   FSH (IU/L) baseline NA 70.1 113.1 3.91 (<0.3–2.4) 8.9 (<0.3–2.4) 3.4 (<0.3–2.4)   FSH (IU/L) GnRH-stimulated NA NA NA 21.62 (1.9–7.6) NA 23.9 (1.9–7.6)   T (ng/dL) baseline 125 (120 ~ 400) 55.0 54.9 <5 (<5) 116 (120–400) <5 (<5)   T (ng/dL) hCG-stimulated 440 (> 200) NA NA 518 (> 200) NA 78 (>200)   DHT (ng/dL) baseline NA NA NA <5 (<5) 20 NA   AMH (ng/mL) baseline NA NA NA 18.3 (43.3–79.3) 19.6 (43.3–79.3) 15.0 (43.3–79.3)  Treatment   Orchidopexy (age) 1 years - - 1 year 7 months Not yet   Genitoplasty (age) 2 years 15 years 15 years 2 years 11 month Not yet   TE (25 mg i.m.) (age) 3 × (6 ~ 10 months) - - No 1 × (5 months)   DHT cream (age) No - - Yes (2 years 8 months) No   PL increment (cm) 2.3→3.9 (+ 1.1 SD) - - Increased (no record) 1.5→3.0 (+1.2 SD)  Central nervous system   Developmental delay No No No Yes No   DQ NA NA NA 63 NA   Seizure No No No No No   ASD/ADHD No No No NA No  Others   Diaphragmatic hernia No No No Yes (L) No NOTE: LH indicates luteinizing hormone; i.m., intramuscular injection; FSH, follicle-stimulating hormone; DHT, dihydrotestosterone; TE, testosterone enanthate; PL, penile length; DQ, development quotient measured with Kyoto Scale of Psychological Development; ASD, autism spectrum disorder; ADHD, attention-deficit hyperactivity disorder. View Large We next investigated the phenotype of the DSD cases with de novo MYRF variants in the four families (Table 1). The three XY patients, Patients 1, 4 and 5, had micro-penis, hypospadias, small testis and cryptorchidism (Supplementary Material, Fig. S2). Laboratory data showed low levels of testosterone (T) and anti-Müllerian hormone (AMH), as well as hypergonadotropism. Intramuscular injection of testosterone enanthate effectively increased the penile length of Patient 4. The low level of AMH is indicative of defects in Sertoli cells. The other findings can be caused by a low level of T during the foetal period, which is cooperatively produced by Sertoli and Leydig cells. One patient showed congenital diaphragmatic hernia (CDH), a cause of which is a defect in the development of pleuroperitoneal folds (13). In the monozygotic twin XX patients of Family 2, uterus, fallopian tube and ovary were hypoplastic, and laboratory data showed hypergonadotropism. These phenotypes can be explained by the defective development of Müllerian derivatives and ovaries. Therefore, the majority of the symptoms observed in the DSD cases of MYRF pathogenic variants can be attributable to defects in the following cells and tissues: Müllerian derivatives, ovaries, Sertoli cells, Leydig cells and pleuroperitoneal folds. Interestingly, all of these cells/tissues are CEDC or partly consist of CEDC (Fig. 2A) and developed via processes of proliferaiton, migration and epithelial-to-mesenchymal transition (EMT) and mesenchymal-to-epithelial transition (MET) (14–18). Thus, the phenotype of the patients suggested that MYRF might play a critical role in the development of CEDC, such as their proliferation, migration, EMT and/or MET. We next investigated the expression patterns of MYRF in foetal gonads. To this end, we analysed public single-cell RNA sequencing (scRNA-seq) data of the human male and female aorta–gonad–mesonephros (AGM) region (4 to 5 wpf) and foetal gonads (7 to 26 wpf) (19). By analysing female gonad data by t-Distributed Stochastic Neighbour Embedding (t-SNE), we identified 11 clusters of possibly different cell identities (C1–11, Fig. 2B). We found that the expression of MYRF was high in C9, 10 and 11 (Fig. 2C and D). We determined the cell identities of each cluster mainly focusing on C9, 10 and 11 as follows (Supplementary Material, Fig. 3). Some cells of C9, C10 and C11 expressed UPK3B, a CE marker (Fig. 2D) (20). C9, C10 and C11 cells also expressed FOXL2, a medullary granulosa cell lineage marker (21), and NR2F2, a putative thecal cell lineage marker (22) (Fig. 2C and D). C9, C10 and C11 rarely expressed LGR5, a cortical granulosa cell lineage marker (22), possibly due to an insufficient sequencing depth (data not shown). Information about CCDC182, another granulosa cell lineage marker (23), was not included in the data, possibly due to the use of old gene annotations. Thus, C9, C10 and C11 could be mixtures of CE, cells of the granulosa cell lineage and cells of the thecal cell lineage. Upon the temporal resolution of each cluster, C9, C10 and C11 mainly consisted of cells at 5–10, 14 and 18–26 wpf and were considered as early, middle and late cells, respectively (Supplementary Material, Fig. S3A). A few cells in the upper part of C9 showed PECAM1 (24) and CDH5 (25) positivity and were considered as endothelial cells (Supplementary Material, Fig. S3M and N). Taking these results together, the properties of cells with high MYRF expression were compatible with CE and CEDC in females, that is, cells of granulosa cell and putative theca cell lineage (Fig. 2A). Figure 2 View largeDownload slide High expression of MYRF in CE, granulosa cell lineage and/or putative theca cell lineage in female gonads. (A) Schematic representation of the relationships among CE and CEDC in foetal gonads. The origin of Leydig cell lineage is partly CE while that of theca cell lineage is not known. Previously reported gene markers for each cell type are as follows: UPK3B (20), FOXL2 (21), NR2F2 (22,27), LGR5 (22), CCDC182 (23), WT1 (26), AMH (26), ARX (28) and CYP17A1 (55). (B–D) scRNA-seq analysis of human female foetal gonads. (B) t-SNE plots showing identities of each cell cluster. Each dot indicates each cell. Dots are coloured based on cell identities determined by t-SNE analysis. Cell type of each cluster was determined by gene markers (C and D) and wpf (Supplementary Material, Fig. S3A). (C) t-SNE plots showing log-normalized expression level of MYRF and marker genes as shown on the right side. Analysed gene is shown above each plot. (D) Bar plot of log-normalized expression level of MYRF and marker genes as shown above each plot. Each bar indicates each cell. Clusters are shown in different colours, so that each cluster is differentiated in the graph. Figure 2 View largeDownload slide High expression of MYRF in CE, granulosa cell lineage and/or putative theca cell lineage in female gonads. (A) Schematic representation of the relationships among CE and CEDC in foetal gonads. The origin of Leydig cell lineage is partly CE while that of theca cell lineage is not known. Previously reported gene markers for each cell type are as follows: UPK3B (20), FOXL2 (21), NR2F2 (22,27), LGR5 (22), CCDC182 (23), WT1 (26), AMH (26), ARX (28) and CYP17A1 (55). (B–D) scRNA-seq analysis of human female foetal gonads. (B) t-SNE plots showing identities of each cell cluster. Each dot indicates each cell. Dots are coloured based on cell identities determined by t-SNE analysis. Cell type of each cluster was determined by gene markers (C and D) and wpf (Supplementary Material, Fig. S3A). (C) t-SNE plots showing log-normalized expression level of MYRF and marker genes as shown on the right side. Analysed gene is shown above each plot. (D) Bar plot of log-normalized expression level of MYRF and marker genes as shown above each plot. Each bar indicates each cell. Clusters are shown in different colours, so that each cluster is differentiated in the graph. Figure 3 View largeDownload slide High expression of MYRF in Sertoli and Leydig cell lineages in male gonads. scRNA-seq analysis of human male foetal gonads. (A) t-SNE plots showing identities of each cell cluster. Each dot indicates each cell. Dots are coloured based on cell identities determined by t-SNE analysis. Cell type of each cluster was determined by gene markers (B, C) and wpf (Supplementary Material, Fig. S4A). (B) t-SNE plots showing log-normalized expression level of MYRF and marker genes as shown on the right side. Analysed gene is shown above each splot. (C) Bar plot of log-normalized expression level of MYRF and marker genes as shown above each plot. Each bar indicates each cell. Clusters are shown in different colours so that each cluster is differentiated in the graph. In addition, early (red) and late (yellow green) Sertoli cells in C6; early (green), middle (blue) and late (purple) Leydig cells in C7; and migrating (red) and other FGC (yellow green) in C1 are coloured differently within each cluster. Figure 3 View largeDownload slide High expression of MYRF in Sertoli and Leydig cell lineages in male gonads. scRNA-seq analysis of human male foetal gonads. (A) t-SNE plots showing identities of each cell cluster. Each dot indicates each cell. Dots are coloured based on cell identities determined by t-SNE analysis. Cell type of each cluster was determined by gene markers (B, C) and wpf (Supplementary Material, Fig. S4A). (B) t-SNE plots showing log-normalized expression level of MYRF and marker genes as shown on the right side. Analysed gene is shown above each splot. (C) Bar plot of log-normalized expression level of MYRF and marker genes as shown above each plot. Each bar indicates each cell. Clusters are shown in different colours so that each cluster is differentiated in the graph. In addition, early (red) and late (yellow green) Sertoli cells in C6; early (green), middle (blue) and late (purple) Leydig cells in C7; and migrating (red) and other FGC (yellow green) in C1 are coloured differently within each cluster. We repeated the above analysis for public data of scRNA-seq on human male foetal gonad. We performed t-SNE analysis and identified eight clusters of possibly different cell identities (C1–8, Fig. 3A). We found that the expression of MYRF was high in C6 diffusely and in some cells of C7 (Fig. 3B). We determined the cell identities of each cluster, mainly focusing on C6 and C7, as follows (Fig. 3B and C and Supplementary Material, Fig. S4) (19). All clusters expressed UPK3B similarly and at a low level (data not shown), possibly indicating that the data contained no CE. C6 expressed WT1 (26) and AMH (26) corresponding to the Sertoli cell lineage (Fig. 3B and C). C7 expressed NR2F2 (27) and ARX (28) and corresponded to the Leydig cell lineage (Fig. 3B and C). Because a sub-cluster of the right part of C6 was mainly from before 10 wpf (Supplementary Material, Fig. S4A), this sub-cluster was considered to consist early cells of the Sertoli cell lineage, while the other part was probably late cells of this lineage. Because a sub-cluster of the upper part of C7 was mainly from before 12 wpf (Supplementary Material, Fig. S4A), the sub-cluster was considered as early cells of the Leydig cell lineage. The few cells at the left part showed positivity for CYP17A1 (29) (Fig. 3B and C), a marker for late cells of the Leydig cell lineage. Thus, the properties of cells with high MYRF expression were compatible with CEDC in males, namely, cells of the Sertoli and Leydig cell lineages (Fig. 2A). Figure 4 View largeDownload slide GO enrichment analysis of putative target genes of Myrf. GO enrichment analysis of putative Myrf target genes in rat differentiating oligodendrocyte progenitor cells. (A) Bar plot of statistical significance (PBH) for enrichment of GO terms of interest related to proliferation, migration and EMT/MET. Bars are coloured according to the categories of GO terms as follows: white for proliferation, grey for migration and black for EMT/MET. The dotted horizontal line indicates PBH of 0.05. (B) Functions of clusters of GO terms with a raw P-value of <0.05 with overlapping genes. Node (white to red circle), GO; node size, correlated with the number of contained genes; node colour, correlated with raw P-value shown as a colour bar at the bottom left; edge (blue lines), overlap of contained genes between GO terms; black oval, GO terms that were similar to each other. Figure 4 View largeDownload slide GO enrichment analysis of putative target genes of Myrf. GO enrichment analysis of putative Myrf target genes in rat differentiating oligodendrocyte progenitor cells. (A) Bar plot of statistical significance (PBH) for enrichment of GO terms of interest related to proliferation, migration and EMT/MET. Bars are coloured according to the categories of GO terms as follows: white for proliferation, grey for migration and black for EMT/MET. The dotted horizontal line indicates PBH of 0.05. (B) Functions of clusters of GO terms with a raw P-value of <0.05 with overlapping genes. Node (white to red circle), GO; node size, correlated with the number of contained genes; node colour, correlated with raw P-value shown as a colour bar at the bottom left; edge (blue lines), overlap of contained genes between GO terms; black oval, GO terms that were similar to each other. We also analysed Myrf expression in another dataset, the Mouse Cell Atlas (30). This dataset contains scRNA-seq data of male and female mouse gonads at E14.5, and we analysed them as we did for the above human foetal gonad data (Supplementary Material, Fig. S5 and S6). In the male and female mouse gonad data, Myrf was rarely expressed, possibly due to a low sequencing depth (Supplementary Material, Fig. S5B and S6B). The cluster in which the mean Myrf expression was the highest was CE (C9 in Supplementary Material, Figs S5 and S6, C8) in both male and female foetal gonads (Supplementary Material, Figs S5C, D and 6C, D). Other clusters expressing Myrf were cells of the granulosa cell lineage and putative thecal cell lineage and FGC in female foetal gonad, while they were early cells of the Leydig cell lineage and cells of the Sertoli cell lineage in male foetal gonad (Supplementary Material, Figs. S5C, D and 6C, D). Myelin regulatory factor protein encoded by MYRF is known as a transcription factor positively regulating gene expression (31,32). We therefore analysed putative target genes of MYRF by using public data of ChIP-seq for Myrf in rat differentiating oligodendrocyte progenitor cells (31). The experiment was performed in a type of cell that differs from those in which we are interested. However, because target genes of transcription factors largely overlap between different types of cell (33), we could partially reveal the target genes of MYRF in affected cells in our patients. We found that Myrf peaks identified by ChIP-seq are highly condensed in regions proximal to transcription start site (TSS; highest within ±10 kb; Supplementary Material, Fig. S7). The results were consistent with a previously reported finding (31) and suggest that these proximal genes are putative targets of Myrf. Thus, we performed a Gene Ontology (GO) enrichment analysis of putative Myrf-target genes whose TSS is close to Myrf peaks (<10 kb) using CHIP-ENRICH (number of putative target genes = 478). As described above, the available ChIP-seq data are derived from oligodendrocyte progenitor cells, and it is expected the target genes are enriched for oligodendritic and other neural processes. Therefore, we first selected a total of 31 GO terms of interest related to proliferation, migration and EMT/MET prior to the analysis (see Materials and Methods) and performed a CHIP-ENRICH analysis. We found that eight and one GO terms for migration and proliferation were significantly [Benjamini–Hochberg-adjusted P-value (PBH) < 0.05] enriched among the putative Myrf target genes, respectively (Fig. 4A). Next, we performed an exploratory enrichment analysis of all GO terms. We identified only three terms with PBH < 0.05: negative regulation of cellular component organization (GO: 0051129), negative regulation of microtubule polymerization or depolymerization (GO: 0031111) and regulation of receptor binding (GO: 1900120) (Supplementary Material, Table S3). To obtain broader insights into the functions of putative Myrf-target genes, we applied a more permissive threshold, raw P < 0.05, and obtained 309 GO terms with nominal significance (Supplementary Material, Table S3). To summarize the functions of these GO terms, we grouped the terms into clusters reflecting their overlapping genes. The two largest clusters were, as expected, related to brain development and development of the nervous system, followed by the regulation of cell migration, regulation of cell-matrix adhesion, regulation of microtubule organization, regulation of response to reactive oxygen species (especially cell death), regulation of JAK-STAT cascade, regulation of kinase (especially MAPK cascade), development of blood vessels and embryonic development (Fig. 4B). Discussion In this study, we identified MYRF as a strong candidate for a gene causative of DSD, for which robust statistically significant results were obtained. The pathogenicity of MYRF damaging variants was supported by multiple lines of evidence: (1) exome-wide significant enrichment of rare LGD variants at MRYF in DSD cases, (2) exome-wide significant enrichment of de novo truncating variants at MYRF in DSD cases, (3) high and selective expression of MYRF in cells and tissues affected in the DSD cases with MYRF pathogenic variants and (4) enrichment of genes related to proliferation and migration among putative target genes of Myrf. Moreover, Pinz et al. (34) recently reported two male cases of DSD and CDH with MYRF truncating variants. While their study alone could not provide robust statistical evidence due to the small number of identified mutation carriers and the lack of a comparison with a control dataset, together with our exome-wide significant findings, there is convincing evidence that MYRF haploinsufficiency causes a type of DSD. Our patients with MYRF variants showed not only gonadal hypoplasia in both 46,XY and 46,XX but also Müllerian duct hypoplasia in 46,XX subjects. The developments of indifferent gonad and Müllerian duct are different biological processes, and gonadal agenesis usually permits Müllerian development. However, agenesis of both gonad and Müllerian duct has been reported in patients with 46,XY (35), 46,XX (36) and 46,X,i(Xq) (37) (since such 46,XY patients with agenesis of both gonad and Müllerian duct exhibit complete female genitalia, it is unlikely that testis tissue is once formed to produce AMH and then regresses towards gonadal hypoplasia). These cases in our study could support the presence of a gene involved in development of both indifferent gonad and Müllerian duct, and the current cases in our study suggest that MYRF is indeed such a candidate. According to the results of our bioinformatic analyses utilizing scRNA-seq and Myrf ChIP-seq data, we found that MYRF/Myrf are highly expressed in CE/CEDC in both human and mouse foetal gonads. GO enrichment analysis of putative Myrf target genes revealed enrichment of GO terms related to migration and proliferation among them. Interpreting these results in an integrated manner, it can be speculated that loss of function of a copy of the MYRF transcription factor leads to the transcriptional dysregulation of genes related to migration and proliferation in CE/CEDC, and then causes the observed phenotypes that can be mostly explained by defects in CE/CEDC development. Nevertheless, this hypothesis could be too simple, and further details of pathomechanisms should be investigated by other approaches, including an analysis of living animals such as heterozygous Myrf knockout mice. We also need to take into account that we used Myrf ChIP-seq data derived from oligodendrocytes (31), which is not a cell type directly affected in our DSD cases. While in this study we report strong evidence indicating a causal relationship between MYRF haploinsufficiency and DSD, MYRF was originally identified as a transcriptional regulator essential for the development and maintenance of central nervous system myelination (32,38). Indeed, MYRF has been reported to be a causative gene of encephalitis/encephalopathy with a reversible splenial lesion (MERS) and MERS-like diseases triggered by febrile illness (39). Considering the expression levels of MYRF in oligodendrocytes and foetal gonads, we found that the Myrf expression level [log10 (mean TPM in cluster +1)] in epithelial cell clusters, which correspond to CE clusters in our analyses, was 4.7 in male foetal gonad and 3.9 in female foetal gonad (data of Mouse Cell Atlas), which was comparable to that in the oligodendrocyte precursor cell cluster (4.8, data of Mouse Cell Atlas). Thus, MYRF can play an important role in CE/CEDC, and it would not be surprising that MYRF damaging variants cause abnormal sex development. Regarding the phenotype–genotype relationship, three out of four DSD cases in our study, and two of two patients in the study of Pinz et al. (34), carry truncating de novo mutations in MYRF. In contrast, the variant reported in families with MERS-like diseases was a missense variant (p.Gln403Arg) that partially decreases the transcriptional activity in vitro (39). Therefore, partial disruption of MYRF function may cause MERS-like diseases that only become apparent after febrile illness and full haploinsufficiency of MYRF would cause DSD. In addition, while in our DSD cases abnormalities of the nervous system were not observed, except in Patient 5 with developmental delay, symptoms that would be induced by febrile illness should be carefully followed. In summary, in this study we demonstrate a causal relationship between MYRF truncating variants and DSD, supported by robust statistical evidence. Our findings can contribute to an improved understanding and diagnosis of DSD. In addition, our integrative bioinformatic analyses provide insights into the putative mechanism by which MYRF haploinsufficiency leads to the observed DSD phenotypes. Finally, the results of this study expand the spectrum of diseases caused by MYRF mutations and clearly point out that MYRF plays important roles not only in the differentiation and maintenance of oligodendrocytes but also in sex development. Materials and Methods Study participants and ethics Informed consent was obtained from all of the study participants. When participants were younger than 20 years, we obtained informed consent from their parents. This study was approved by the institutional review board of Yokohama City University School of Medicine (Yokohama, Japan). Trained obstetricians, gynaecologists and/or paediatricians performed the DSD diagnosis and the detailed clinical assessment of the patients. Control samples in the following case–control analysis were from unrelated healthy parents who had children suffering from rare diseases. Whole-exome sequencing WES was performed for 26 DSD and 2625 control samples as previously described (37–40). In brief, genomic DNA was captured with a SureSelect Human All Exon V4, V5 or V6 Kit (Agilent Technologies, Santa Clara, CA, USA): V4 in 8, V5 in 16 and V6 in 2 in 26 DSD samples while V4 in 144, V5 in 1567 and V6 in 914 samples of 2625 control samples. Sequencing was performed on an Illumina HiSeq 2500 (Illumina, San Diego, CA, USA) with 101 bp paired-end reads. Reads were aligned to the human reference genome (GRCh37.1/hg19) using Novoalign v3.02.13. Polymerase chain reaction (PCR) duplicates were removed using Picard. Local realignments around indels and base quality score recalibration were performed with the Genome Analysis Toolkit (GATK) 3.7-0, and analysis-ready binary alignment map (BAM) files were generated (41). From the BAM files, variants in all case and control samples were jointly genotyped with GATK HaplotypeCaller and hard-filtered in accordance with GATK Best Practice (42). To remove possible bias from different proportions of capture kits used in DSD and control groups, we selected variants only at bait regions overlapped among V4, V5 and V6 capture kits with a margin of 100 bp using merge and intersect function of bedtools. Variants were annotated using dbNSFP in ANNOVAR and snpEff (43–45). Case–control analysis of rare variant burden The burden of rare damaging variants at each protein-coding gene was compared between 26 DSD cases and 2625 controls. We defined damaging variants in the following two ways: LGD or LGD + D-mis at each gene. We defined variants of HIGH impact in snpEff (e.g. stop_gained, frameshift_variant and splice_donor/acceptor_variant) as LGD and missense variants with a SIFT (8) score of <0.05 and a PolyPhen2 HVAR (7) score of >0.90 as D-mis. We defined rare variants as follows: (1) minor allele frequency (MAF) of <0.0005 in databases of variants in general populations: ESP6500, gnomAD (46) and ToMMo3.5KJPN (47) and (2) not existing in our intramural databases of 575 Japanese controls. We compared the burden of these variants in each gene between DSD cases and controls by a collapsing method (48) followed by Fisher’s exact test using RVTESTS (49). We conservatively defined the threshold P-value for statistical significance by applying Bonferroni correction with the total number of tests in our LGD only and LGD + D-mis models. After excluding genes with no rare damaging variant among the overall cohort (26 DSD cases and 2625 controls), to which Fisher’s exact test cannot be applied, there were 9217 and 15 177 testable genes in LGD only and LGD + D-mis models, respectively. With these numbers of testable genes, we defined the threshold for exome-wide significance as 2.05E−06 [= 0.05/(9217 + 15 177)]. Analysis of burden of de novo variants The burden of the de novo truncating variant at MYRF in 26 DSD cases was analysed under a null model of mutation (6) that compares observed and expected numbers of de novo variants in a gene of interest. Briefly, numbers of de novo variants are expected to follow a Poisson distribution with λ (expected number of de novo variants) calculated as follows: (rate of de novo truncating variants at gene of interest per allele per generation) × (number of tested alleles). Samocha et al. (6) reported the probabilities of de novo nonsense, frameshift, canonical splice-site and missense variants per generation in a copy of MYRF as 1.12E-06, 2.19E-06, 1.06E-06 and 3.29E-05 and thereby the probability of LGD variants or LGD + missense variants can be considered as the sum of these probabilities (4.36E-06 or 3.73E-05, respectively). Thus, the expected count of de novo MYRF LGD variants in 26 cases or LGD + missense variants in 94 cases is expected to follow a Poisson distribution with a mean of 2.27E-04 or 7.01E-03, respectively (= 4.36E-06 × 26 × 2 or 3.73E-05 × 94 × 2). Because we observed three de novo LGD variants in MYRF in 26 DSD cases and four de novo LGD + missense variants in 94 DSD cases, the probability of three or more counts in the Poisson distribution was calculated with the POISSON function of Microsoft Excel. The threshold for statistical significance was 2.05E-06 as described in the above case–control analysis of the rare variant burden. Sanger sequencing The identified MYRF variants in the probands and the family members were analysed by Sanger sequencing using standard methods. Briefly, PCR products were purified with alkaline phosphatase and exonuclease 1 and sequenced with BigDye Terminator v3.1 Cycle Sequencing kit (Applied Biosystems, Foster City, CA, USA) on a 3500 DNA Sequencing Analyzer (Life Technologies, Carlsbad, CA, USA). Primer sequences are available on request. scRNA-seq analysis We reanalysed two public sets of data of scRNA-seq: (1) data of the human AGM region at 5 weeks post-fertilization (wpf) and foetal gonads from 7 to 26 wpf (19) and (2) data of mouse foetal gonads (30). We analysed the data with Seurat (50) in accordance with the developer’s guide. Briefly, we retained genes expressed in at least three cells and retained cells with 200–2500 detected genes using the FilterCells function. We normalized the gene expression by the total expression in each cell, multiplied it by 10 000 and log10-transformed it (hereafter, log-normalized expression level) using the NormalizeData function. We extracted genes with highly variable expression using the FindVariableGenes function and the following thresholds: mean of 0.1-3 and dispersion of log-normalized expression level of >0.5. We linear-regressed out the log-normalized expression level using the number of unique molecular identifiers in each cell, and scaled and centred residuals (hereafter, scaled expression level). We performed principal component analysis of the scaled expression level of the genes with highly variable expression and identified several PCs in each set of data that had a clearly larger standard deviation than the remaining PCs in PCElbowPlot function. We performed t-SNE using the RunTSNE function and the PCs with a large standard deviation as input and plotted the data for each cell in two dimensions. ChIP-seq analysis We reanalysed public data of ChIP-seq for Myrf (28). Briefly, in the originally reported work, cultured rat oligodendrocyte progenitor cells were transfected with vectors expressing either Myc-MYRF or untagged MYRF. After 48 h in differentiating conditions, ChIP-seq using Myc-tag antibody was performed, and 2102 and 17 peaks were detected by MACS (51) in Myc-MYRF and untagged MYRF experiments, respectively. Because all of the 17 peaks were detected among the 2102 peaks, there were 2085 peaks specific to the Myc-MYRF experiments. Using information on these 2085 peaks, we analysed the enrichment of biological process-related GO terms among the putative target genes of Myrf. For this purpose, we used CHIP-ENRICH (52) with the default settings. The Locus Definition parameter, which defines the threshold of distance between a Myrf peak and a TSS of a putative target gene, was determined as ‘<10 kb’, according to the relationship between the Myrf peaks and TSS. In this analysis, we analysed the enrichment of GO terms related to proliferation (number of terms, 14), migration (number of terms, 12) and EMT/MET (number of terms, 5), which were selected prior to the analysis and corrected the obtained raw P-values using the Benjamini–Hochberg method. To summarize the functions of all GO terms with a raw P-value of <0.05 in the CHIP-ENRICH results, we clustered the GO terms whose genes overlapped with >0.7 overlap coefficient using the EnrichmentMap (53) plugin of CytoScape v3.6.1 (54). To prevent the overclustering of unrelated GO terms, we excluded GO terms that contained >1000 genes. The distance from Myrf peaks to the nearest TSS outputted by the CHIP-ENRICH analysis was used to generate plots of the relationship between TSS and Myrf peaks. Acknowledgements We thank all of the participants for their cooperation in this research. We also thank Ms K. Takabe, Mr T. Miyama, Ms N. Watanabe, Ms M. Sato, Mr S. Nakamura and Ms S. Sugimoto at the Department of Human Genetics, Yokohama City University Graduate School of Medicine, for their technical assistance. We are also grateful to Edanz (www.edanzediting.co.jp) for editing the English text of a draft of this manuscript. Conflict of Interest statement. None declared. Funding Japan Agency for Medical Research and Development (AMED) (grant numbers JP18ek0109280, JP18dm0107090, JP18ek0109301, JP18ek0109348 and JP18kk020500 to N. Matsumoto); Japan Society for the Promotion of Science (JSPS) KAKENHI (grant numbers JP17H01539 to N. Matsumoto, JP16H05160 to H.S., JP16H05357 to N. Miyake, JP16H06254 to A.T., JP15K10367 to M.N., JP17K10080 to S.M., JP17K15630 to T.M., JP17H06994 to A.F.; Ministry of Health, Labour and Welfare (to N. Matsumoto); the Takeda Science Foundation (to N. Miyake, H.S, N. Matsumoto); and The Ichiro Kanehara Foundation for the Promotion of Medical Science & Medical Care (S.M.). References 1 Lee , P.A. , Houk , C.P. , Ahmed , S.F. and Hughes , I.A. ( 2006 ) Consensus statement on management of intersex disorders. International Consensus Conference on Intersex . Pediatrics , 118 , e488 – e500 . Google Scholar Crossref Search ADS PubMed WorldCat 2 Bashamboo , A. and McElreavey , K. ( 2013 ) Gene mutations associated with anomalies of human gonad formation . Sex. Dev. , 7 , 126 – 146 . Google Scholar Crossref Search ADS PubMed WorldCat 3 Dong , Y. , Yi , Y. , Yao , H. , Yang , Z. , Hu , H. , Liu , J. , Gao , C. , Zhang , M. , Zhou , L. , Asan , L. , Yi , X. and Liang , Z. ( 2016 ) Targeted next-generation sequencing identification of mutations in patients with disorders of sex development . BMC Med. Genet. , 17 , 23 . Google Scholar Crossref Search ADS PubMed WorldCat 4 Baxter , R.M. , Arboleda , V.A. , Lee , H. , Barseghyan , H. , Adam , M.P. , Fechner , P.Y. , Bargman , R. , Keegan , C. , Travers , S. , Schelley , S. et al. ( 2015 ) Exome sequencing for the diagnosis of 46,XY disorders of sex development . J. Clin. Endocrinol. Metab. , 100 , E333 – E344 . Google Scholar Crossref Search ADS PubMed WorldCat 5 Richards , S. , Aziz , N. , Bale , S. , Bick , D. , Das , S. , Gastier-Foster , J. , Grody , W.W. , Hegde , M. , Lyon , E. , Spector , E. , Voelkerding , K. and Rehm , H.L. ( 2015 ) Standards and guidelines for the interpretation of sequence variants: a joint consensus recommendation of the American College of Medical Genetics and Genomics and the Association for Molecular Pathology . Genet. Med. , 17 , 405 – 424 . Google Scholar Crossref Search ADS PubMed WorldCat 6 Samocha , K.E. , Robinson , E.B. , Sanders , S.J. , Stevens , C. , Sabo , A. , McGrath , L.M. , Kosmicki , J.A. , Rehnstrom , K. , Mallick , S. , Kirby , A. et al. ( 2014 ) A framework for the interpretation of de novo mutation in human disease . Nat. Genet. , 46 , 944 – 950 . Google Scholar Crossref Search ADS PubMed WorldCat 7 Adzhubei , I.A. ( 2010 ) A method and server for predicting damaging missense mutations . Nat. Methods , 7 , 248 – 249 . Google Scholar Crossref Search ADS PubMed WorldCat 8 Ng , P.C. and Henikoff , S. ( 2003 ) SIFT: predicting amino acid changes that affect protein function . Nucleic Acids Res. , 31 , 3812 – 3814 . Google Scholar Crossref Search ADS PubMed WorldCat 9 Schwarz , J.M. , Cooper , D.N. , Schuelke , M. and Seelow , D. ( 2014 ) MutationTaster2: mutation prediction for the deep-sequencing age . Nat. Methods , 11 , 361 – 362 . Google Scholar Crossref Search ADS PubMed WorldCat 10 Kircher , M. , Witten , D.M. , Jain , P. , O'Roak , B.J. , Cooper , G.M. and Shendure , J. ( 2014 ) A general framework for estimating the relative pathogenicity of human genetic variants . Nat. Genet. , 46 , 310 – 315 . Google Scholar Crossref Search ADS PubMed WorldCat 11 Zhen , X. , Li , B. , Hu , F. , Yan , S. , Meloni , G. , Li , H. and Shi , N. ( 2017 ) Crystal structure of the DNA-binding domain of myelin-gene regulatory factor . Sci. Rep. , 7 , 3696 . Google Scholar Crossref Search ADS PubMed WorldCat 12 Firth , H.V. , Richards , S.M. , Bevan , A.P. , Clayton , S. , Corpas , M. , Rajan , D. , Van Vooren , S. , Moreau , Y. , Pettett , R.M. and Carter , N.P. ( 2009 ) DECIPHER: database of chromosomal imbalance and phenotype in humans using Ensembl resources . Am. J. Hum. Genet. , 84 , 524 – 533 . Google Scholar Crossref Search ADS PubMed WorldCat 13 Kardon , G. , Ackerman , K.G. , McCulley , D.J. , Shen , Y. , Wynn , J. , Shang , L. , Bogenschutz , E. , Sun , X. and Chung , W.K. ( 2017 ) Congenital diaphragmatic hernias: from genes to mechanisms to therapies . Dis. Model. Mech. , 10 , 955 – 970 . Google Scholar Crossref Search ADS PubMed WorldCat 14 Orvis , G.D. and Behringer , R.R. ( 2007 ) Cellular mechanisms of Mullerian duct formation in the mouse . Dev. Biol. , 306 , 493 – 504 . Google Scholar Crossref Search ADS PubMed WorldCat 15 Guioli , S. , Sekido , R. and Lovell-Badge , R. ( 2007 ) The origin of the Mullerian duct in chick and mouse . Dev. Biol. , 302 , 389 – 398 . Google Scholar Crossref Search ADS PubMed WorldCat 16 Carmona , R. , Canete , A. , Cano , E. , Ariza , L. , Rojas , A. and Munoz-Chapuli , R. ( 2016 ) Conditional deletion of WT1 in the septum transversum mesenchyme causes congenital diaphragmatic hernia in mice . Elife , 5 . doi: https://doi.org/10.7554/eLife.16009 . WorldCat 17 Piprek , R.P. ( 2016 ) Molecular Mechanisms of Cell Differentiation in Gonad Development . Springer Publishing , NY . 18 Ariza , L. , Carmona , R. , Canete , A. , Cano , E. and Munoz-Chapuli , R. ( 2016 ) Coelomic epithelium-derived cells in visceral morphogenesis . Dev. Dyn. , 245 , 307 – 322 . Google Scholar Crossref Search ADS PubMed WorldCat 19 Li , L. , Dong , J. , Yan , L. , Yong , J. , Liu , X. , Hu , Y. , Fan , X. , Wu , X. , Guo , H. , Wang , X. et al. ( 2017 ) Single-cell RNA-seq analysis maps development of human germline cells and gonadal niche interactions . Cell Stem Cell , 20 , 858 – 873 e854 . Google Scholar Crossref Search ADS PubMed WorldCat 20 Rudat , C. , Grieskamp , T. , Rohr , C. , Airik , R. , Wrede , C. , Hegermann , J. , Herrmann , B.G. , Schuster-Gossler , K. and Kispert , A. ( 2014 ) Upk3b is dispensable for development and integrity of urothelium and mesothelium . PLoS One , 9 , e112112. doi: https://doi.org/10.1371/journal.pone.0112112 . WorldCat 21 Mork , L. , Maatouk , D.M. , McMahon , J.A. , Guo , J.J. , Zhang , P. , McMahon , A.P. and Capel , B. ( 2012 ) Temporal differences in granulosa cell specification in the ovary reflect distinct follicle fates in mice . Biol. Reprod. , 86 , 37 . Google Scholar Crossref Search ADS PubMed WorldCat 22 Rastetter , R.H. , Bernard , P. , Palmer , J.S. , Chassot , A.A. , Chen , H. , Western , P.S. , Ramsay , R.G. , Chaboissier , M.C. and Wilhelm , D. ( 2014 ) Marker genes identify three somatic cell types in the fetal mouse ovary . Dev. Biol. , 394 , 242 – 252 . Google Scholar Crossref Search ADS PubMed WorldCat 23 Lee , H.J. , Pazin , D.E. , Kahlon , R.S. , Correa , S.M. and Albrecht , K.H. ( 2009 ) Novel markers of early ovarian pre-granulosa cells are expressed in an Sry-like pattern . Dev. Dyn. , 238 , 812 – 825 . Google Scholar Crossref Search ADS PubMed WorldCat 24 Stockinger , H. , Gadd , S.J. , Eher , R. , Majdic , O. , Schreiber , W. , Kasinrerk , W. , Strass , B. , Schnabl , E. and Knapp , W. ( 1990 ) Molecular characterization and functional analysis of the leukocyte surface protein CD31 . J. Immunol. , 145 , 3889 – 3897 . Google Scholar PubMed WorldCat 25 Gory , S. , Vernet , M. , Laurent , M. , Dejana , E. , Dalmon , J. and Huber , P. ( 1999 ) The vascular endothelial-cadherin promoter directs endothelial-specific expression in transgenic mice . Blood , 93 , 184 – 192 . Google Scholar PubMed WorldCat 26 de Santa Barbara , P. , Moniot , B. , Poulat , F. and Berta , P. ( 2000 ) Expression and subcellular localization of SF-1, SOX9, WT1, and AMH proteins during early human testicular development . Dev. Dyn. , 217 , 293 – 298 . Google Scholar Crossref Search ADS PubMed WorldCat 27 van den Driesche , S. , Walker , M. , McKinnell , C. , Scott , H.M. , Eddie , S.L. , Mitchell , R.T. , Seckl , J.R. , Drake , A.J. , Smith , L.B. , Anderson , R.A. and Sharpe , R.M. ( 2012 ) Proposed role for COUP-TFII in regulating fetal Leydig cell steroidogenesis, perturbation of which leads to masculinization disorders in rodents . PLoS One , 7 , e37064. doi: https://doi.org/10.1371/journal.pone.0037064 . WorldCat 28 Miyabayashi , K. , Katoh-Fukui , Y. , Ogawa , H. , Baba , T. , Shima , Y. , Sugiyama , N. , Kitamura , K. and Morohashi , K. ( 2013 ) Aristaless related homeobox gene, Arx, is implicated in mouse fetal Leydig cell differentiation possibly through expressing in the progenitor cells . PLoS One , 8 , e68050. doi: https://doi.org/10.1371/journal.pone.0068050 . WorldCat 29 Miyabayashi , K. , Shima , Y. , Inoue , M. , Sato , T. , Baba , T. , Ohkawa , Y. , Suyama , M. and Morohashi , K.I. ( 2017 ) Alterations in fetal Leydig cell gene expression during fetal and adult development . Sex. Dev. , 11 , 53 – 63 . Google Scholar Crossref Search ADS PubMed WorldCat 30 Han , X. , Wang , R. , Zhou , Y. , Fei , L. , Sun , H. , Lai , S. , Saadatpour , A. , Zhou , Z. , Chen , H. , Ye , F. et al. ( 2018 ) Mapping the mouse cell atlas by microwell-seq . Cell , 172 , 1091 – 1107 . Google Scholar Crossref Search ADS PubMed WorldCat 31 Bujalka , H. , Koenning , M. , Jackson , S. , Perreau , V.M. , Pope , B. , Hay , C.M. , Mitew , S. , Hill , A.F. , Lu , Q.R. , Wegner , M. et al. ( 2013 ) MYRF is a membrane-associated transcription factor that autoproteolytically cleaves to directly activate myelin genes . PLoS Biol , 11 , e1001625. doi: https://doi.org/10.1371/journal.pbio.1001625 . WorldCat 32 Emery , B. , Agalliu , D. , Cahoy , J.D. , Watkins , T.A. , Dugas , J.C. , Mulinyawe , S.B. , Ibrahim , A. , Ligon , K.L. , Rowitch , D.H. and Barres , B.A. ( 2009 ) Myelin gene regulatory factor is a critical transcriptional regulator required for CNS myelination . Cell , 138 , 172 – 185 . Google Scholar Crossref Search ADS PubMed WorldCat 33 Handstad , T. , Rye , M. , Mocnik , R. , Drablos , F. and Saetrom , P. ( 2012 ) Cell-type specificity of ChIP-predicted transcription factor binding sites . BMC Genomics , 13 , 372 . Google Scholar Crossref Search ADS PubMed WorldCat 34 Pinz , H. , Pyle , L.C. , Li , D. , Izumi , K. , Skraban , C. , Tarpinian , J. , Braddock , S.R. , Telegrafi , A. , Monaghan , K.G. , Zackai , E. and Bhoj , E.J. ( 2018 ) De novo variants in myelin regulatory factor (MYRF) as candidates of a new syndrome of cardiac and urogenital anomalies . Am. J. Med. Genet. A , 176 , 969 – 972 . Google Scholar Crossref Search ADS PubMed WorldCat 35 Sarto , G.E. and Opitz , J.M. ( 1973 ) The XY gonadal agenesis syndrome . J. Med. Genet. , 10 , 288 – 293 . Google Scholar Crossref Search ADS PubMed WorldCat 36 Levinson , G. , Zarate , A. , Guzman-Toledano , R. , Canales , E.S. and Jimenez , M. ( 1976 ) An XX female with sexual infantilism, absent gonads, and lack of Mullerian ducts . J. Med. Genet. , 13 , 68 – 69 . Google Scholar Crossref Search ADS PubMed WorldCat 37 De Leon , F.D. , Hersh , J.H. , Sanfilippo , J.S. , Schikler , K.N. and Yen , F.F. ( 1984 ) Gonadal and Mullerian duct agenesis in a girl with 46,X,i(Xq) . Obstet. Gynecol. , 63 , 81s – 83s . Google Scholar PubMed WorldCat 38 Koenning , M. , Jackson , S. , Hay , C.M. , Faux , C. , Kilpatrick , T.J. , Willingham , M. and Emery , B. ( 2012 ) Myelin gene regulatory factor is required for maintenance of myelin and mature oligodendrocyte identity in the adult CNS . J. Neurosci. , 32 , 12528 – 12542 . Google Scholar Crossref Search ADS PubMed WorldCat 39 Kurahashi , H. , Azuma , Y. , Masuda , A. , Okuno , T. , Nakahara , E. , Imamura , T. , Saitoh , M. , Mizuguchi , M. , Shimizu , T. , Ohno , K. and Okumura , A. ( 2018 ) MYRF is associated with encephalopathy with reversible myelin vacuolization . Ann. Neurol. , 83 , 98 – 106 . Google Scholar Crossref Search ADS PubMed WorldCat 40 Hamanaka , K. , Miyatake , S. , Zerem , A. , Lev , D. , Blumkin , L. , Yokochi , K. , Fujita , A. , Imagawa , E. , Iwama , K. , Nakashima , M. et al. ( 2018 ) Expanding the phenotype of IBA57 mutations: related leukodystrophy can remain asymptomatic . J. Hum. Genet. , 63 , 1223 – 1229 . Google Scholar Crossref Search ADS PubMed WorldCat 41 McKenna , A. , Hanna , M. , Banks , E. , Sivachenko , A. , Cibulskis , K. , Kernytsky , A. , Garimella , K. , Altshuler , D. , Gabriel , S. , Daly , M. and DePristo , M.A. ( 2010 ) The genome analysis toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data . Genome Res. , 20 , 1297 – 1303 . Google Scholar Crossref Search ADS PubMed WorldCat 42 DePristo , M.A. , Banks , E. , Poplin , R. , Garimella , K.V. , Maguire , J.R. , Hartl , C. , Philippakis , A.A. , del Angel , G. , Rivas , M.A. , Hanna , M. et al. ( 2011 ) A framework for variation discovery and genotyping using next-generation DNA sequencing data . Nat. Genet. , 43 , 491 – 498 . Google Scholar Crossref Search ADS PubMed WorldCat 43 Wang , K. , Li , M. and Hakonarson , H. ( 2010 ) ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data . Nucleic Acids Res. , 38 , e164 . Google Scholar Crossref Search ADS PubMed WorldCat 44 Cingolani , P. , Platts , A. , Wang le , L. , Coon , M. , Nguyen , T. , Wang , L. , Land , S.J. , Lu , X. and Ruden , D.M. ( 2012 ) A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3 . Fly (Austin) , 6 , 80 – 92 . Google Scholar Crossref Search ADS PubMed WorldCat 45 Liu , X. , Jian , X. and Boerwinkle , E. ( 2011 ) dbNSFP: a lightweight database of human nonsynonymous SNPs and their functional predictions . Hum. Mutat. , 32 , 894 – 899 . Google Scholar Crossref Search ADS PubMed WorldCat 46 Lek , M. ( 2016 ) Analysis of protein-coding genetic variation in 60,706 humans . Nature , 536 , 285 – 291 . Google Scholar Crossref Search ADS PubMed WorldCat 47 Nagasaki , M. , Yasuda , J. , Katsuoka , F. , Nariai , N. , Kojima , K. , Kawai , Y. , Yamaguchi-Kabata , Y. , Yokozawa , J. , Danjoh , I. , Saito , S. et al. ( 2015 ) Rare variant discovery by deep whole-genome sequencing of 1,070 Japanese individuals . Nat. Commun. , 6 , 8018 . Google Scholar Crossref Search ADS PubMed WorldCat 48 Li , B. and Leal , S.M. ( 2008 ) Methods for detecting associations with rare variants for common diseases: application to analysis of sequence data . Am. J. Hum. Genet. , 83 , 311 – 321 . Google Scholar Crossref Search ADS PubMed WorldCat 49 Zhan , X. , Hu , Y. , Li , B. , Abecasis , G.R. and Liu , D.J. ( 2016 ) RVTESTS: an efficient and comprehensive tool for rare variant association analysis using sequence data . Bioinformatics , 32 , 1423 – 1426 . Google Scholar Crossref Search ADS PubMed WorldCat 50 Butler , A. , Hoffman , P. , Smibert , P. , Papalexi , E. and Satija , R. ( 2018 ) Integrating single-cell transcriptomic data across different conditions, technologies, and species . Nat. Biotechnol. , 36 , 411 – 420 . Google Scholar Crossref Search ADS PubMed WorldCat 51 Zhang , Y. , Liu , T. , Meyer , C.A. , Eeckhoute , J. , Johnson , D.S. , Bernstein , B.E. , Nusbaum , C. , Myers , R.M. , Brown , M. , Li , W. and Liu , X.S. ( 2008 ) Model-based analysis of ChIP-seq (MACS) . Genome Biol. , 9 , R137 . Google Scholar Crossref Search ADS PubMed WorldCat 52 Welch , R.P. , Lee , C. , Imbriano , P.M. , Patil , S. , Weymouth , T.E. , Smith , R.A. , Scott , L.J. and Sartor , M.A. ( 2014 ) ChIP-Enrich: gene set enrichment testing for ChIP-seq data . Nucleic Acids Res. , 42 , e105 . Google Scholar Crossref Search ADS PubMed WorldCat 53 Merico , D. , Isserlin , R. , Stueker , O. , Emili , A. and Bader , G.D. ( 2010 ) Enrichment map: a network-based method for gene-set enrichment visualization and interpretation . PLoS One , 5 , e13984. WorldCat 54 Shannon , P. , Markiel , A. , Ozier , O. , Baliga , N.S. , Wang , J.T. , Ramage , D. , Amin , N. , Schwikowski , B. and Ideker , T. ( 2003 ) Cytoscape: a software environment for integrated models of biomolecular interaction networks . Genome Res. , 13 , 2498 – 2504 . Google Scholar Crossref Search ADS PubMed WorldCat © The Author(s) 2019. Published by Oxford University Press. All rights reserved. For Permissions, please email: journals.permissions@oup.com This article is published and distributed under the terms of the Oxford University Press, Standard Journals Publication Model (https://academic.oup.com/journals/pages/open_access/funder_policies/chorus/standard_publication_model)

Journal

Human Molecular Genetics – Oxford University Press

Published: Jul 15, 2019

There are no references for this article.