Exome sequencing and analysis

DNA was isolated from peripheral blood using Gentra Puregene Blood Kit (Qiagen, USA) or saliva collected using ORAgene ORG-500 (DNAgenoteck, Canada). Sequencing libraries and exome capture was done for each sample following manufacturer’s protocols for SureSelect All Exon 50 Mb capture kit (Agilent Technologies) and Nextera Rapid Capture (Illumina, USA). Sequencing was performed on an Illumina HiSeq2500 as 50-100 bp paired-end run at the UCLA Clinical Genomics Center.

The sequence reads, FASTQ files, were aligned to the human reference genome (GRCh37/hg19 Feb. 2009 assembly) using BWA (Burrows-Wheeler Alignment tool) [27] and Novoalign (novocraft.com). The output BAM files were sorted and merged, and PCR duplicates were removed using Picard. INDEL (insertion and deletion) realignment and recalibration was performed using Genome Analysis Tool Kit (GATK) (broadinstitute.org). Both single-nucleotide variants (SNVs) and small INDELs were called within the Ensembl coding exonic intervals ± 2 bp using GATK’s Unified Genotyper, then recalibrated and filtered using GATK variant-quality score recalibration and variant filtration tools. All high-quality variants were annotated using SNP&Variation Suite and VarSeq—variant filtration and annotation software (Golden Helix, USA). All variants were filtered by a minor allele frequency (MAF) of < 1% and intersected with the DSD gene list to identify mutations in known DSD genes. The list is comprised of a primary gene list of well-annotated genes involved in sex determination and differentiation [17], as well as a secondary list of genes that are more loosely associated with sex development, e.g., their OMIM (Online Mendelian Inheritance of Man) description contains sex development keywords.

The variants identified by exome sequencing were classified into causative or likely causative variants following the recommendations of the American College of Medical Genetics and Genomics [28]. All other variants with minor allele frequency below 1% were classified as variants of unknown significance (VUS). To assess previously unreported missense variants, we used two in silico algorithms SIFT [29] and PolyPhen [30] to predict the pathogenicity of a missense variant based on conservation of the amino acid across species, the physical characteristics of the altered amino acid, and the possible impact on protein structure and function. All variants with low quality scores were validated by Sanger sequencing [31].

Animal care and dissections

The C57BL/6J and C57BL/6J-Y ^POS animals were housed at the UCLA Animal Care Facility following the guidelines of the University of California, Los Angeles, Division of Laboratory Animal Medicine. All experiments were approved by the Institutional Animal Care and Research Committees of UCLA. Wild-type C57BL/6J males and females used for breeding were purchased from the Jackson Laboratory (Bar Harbor, USA), which is fully accredited by the American Association for Accreditation of Laboratory Animal Care.

We have previously identified a 1.5-Mb congenic region on chromosome 11 that confers 80% protection from B6-Y ^POS sex reversal in the heterozygous state (B6-110h-Y ^POS ) and complete protection in the homozygous state (B6-110H-Y ^POS ) [25]. This protective region allows for continual maintenance of subfertile poschiavinus male mice as a breeding colony, with an option of generating unprotected B6-Y ^POS males by mating heterozygous B6-110h-Y ^POS males with wild-type (WT) B6 females. Overnight mating was performed using either the wild-type (WT) B6 or B6-110h-Y ^POS (protected from sex reversal) males and WT B6 females. Dissections were performed at E11.5; the gonads were separated from the mesonephros and placed in RNA stabilizing solution RNAlater (Ambion). DNA was extracted from the rest of the embryos for genotyping. Chromosomal sex was determined using a single primer pair for X-linked Smc-x gene (330 bp) and the Y-linked Smc-y gene (301 bp) (forward: 5′CCGCTGCCAAATTCTTTGG3′; reverse: 5′TGAAGCTTTTGGCTTTGAG3′). The presence of the Y ^POS chromosome was determined by a SNP between Y^B6 and Y ^POS Sry gene using the primer sets 5′TGAATGCATTTATGGTGTGGTC3′; 5′AGCTTTGCTGGTTTTTGGAGTA3′. Immomix Red (Bioline, UK). Presence or absence of the 110 h protective region in B6-Y ^POS males was checked by Sanger sequencing of two regions 11-10 and 11-11 containing SNP rs27019103 (5′AAAGTGTGCTTCCCAGGAGA3′; 5′CCTCTCCCTCAACCCCTAAG3′) and SNP rs28240850 (5′CCACAGCTGGAGGTAGGGTA3′; 5′CCTAAGATGCCATGGGAAGA3′) respectively [25].

Total RNA was isolated from combined embryonic gonadal tissue (50–70 gonads per group) using Qiagen RNeasy Kit (Qiagen) following manufacturer’s guidelines. RNA quality was assessed by Agilent 2100 Bioanalyzer (Agilent Technologies). All samples were required to have RNA integrity scores (RIN) greater than 8.

Experiments on AmhCre Sox9floxflox mice were carried out in strict adherence with the recommendations in the Australian code of practice for the care and use of animals for scientific purposes from the National Health and Medical Research Council.

E13.5 gonads were separated from the mesonephros and total RNA was isolated using the RNeasy Kit (Qiagen), as described in Rahmoun et al. [32]. Embryos were genotyped, and the RNA from the six gonads was pooled into wild-type XY, XX, or XY AmhCre Sox9floxflox (Sox9 knockout). This was repeated in three biological replicates; protocols are detailed in Rahmoun et al. [32].

RNA sequencing and expression analysis

RNA from each sample was submitted to the UCLA Neuroscience Genomic Core (UNGC) for library preparation and sequencing. Library preparation was performed using TruSeq Stranded Total RNA kit (Illumina) with Poly-A selection following manufacturer’s guidelines. Sequencing was performed on HiSeq 2500 (Illumina) with 69 bp paired-end run on a rapid flowcell capable of generating 150 M reads per lane. Four samples were multiplexed and sequenced over two rapid lanes with each sample receiving approximately 75 million reads with > 85% map rate.

The generated sequencing reads were aligned to the mouse genome, version mm10 with STAR [33]. Transcript abundance was assessed by Cufflinks (v2.1.1) [34], using a GTF file based on Ensembl mouse NCBI37. Differential expression analysis was based on fold change differences greater than 1.5 between the groups being compared. Differentially expressed genes were split into two categories: underexpressed and overexpressed in B6-Y ^POS males. Both categories were separately subjected to pathway enrichment analysis using Gene Ontology Consortium [35].

To analyze the RNA from AmhCre Sox9floxflox XY gonads and wild-type XY and XX gonads, libraries were generated using the NuGEN Mondrian Technology and SPIA amplification methodology, and the data was processed and aligned to the mouse genome (Ensembl version 38.77) as described by Rahmoun et al. [32]. To eliminate composition biases, the trimmed mean of M values (TMM) method was used for normalization between the samples [36]. The adjusted P value of 0.05 was used to assess which genes were differentially expressed between XY and AmhCre Sox9floxflox XY (Sox9 KO). Graphs of gene expression were made using GraphPad Prism.

Quantitative PCR

Reverse transcription of RNA to cDNA was performed using Tetro cDNA Synthesis Kit (Bioline, UK) following manufacturer’s protocol. The primer sequences used are detailed in Additional file 1: Table S1. Primers were designed using autoprime software (autoprime.de) and spanned exon-exon junctions for optimal RNA quantification. cDNA was quantified using QuBit HS (Invitrogen) for double-stranded DNA, and a total of 3 ng of cDNA was used per sample for amplification. qPCR was carried out in duplicates using SensiFAST™ SYBR No-ROX Kit (Bioline, UK) by DNA Engine Opticon® 2 real-time PCR detection system (BioRad, USA). Reaction conditions were as follows: 95 °C for 10 min, then 40 cycles of 95 °C for 15 s, 60–64 °C (see Additional file 1: Table S1) for 10 s, and 72 °C for 15 s. Data was analyzed via Opticon Monitor Software (BioRad). Standard curves were generated from a mix of cDNA of all tested samples with five iterations of 1:4 dilutions. Average cycle threshold values (Ct) for each gene/sample were determined based on two replicates. Complementary DNA amounts were estimated based on Ct values and linear equation y = mx + b (where y is the Ct value, m is the slope, x is the cDNA amount, and b is the intercept).

Immunohistochemistry

Fbln2 expression in the embryonic gonads at E12.5 was assessed using immunohistochemistry following the experimental design of Wilhelm et al. [37], using the anti-Fbln2 rabbit polyclonal, sc-30176 (Santa Cruz Biotechnology) antibody. Topro (Invitrogen) was used to counterstain nuclei. All images were taken on a Zeiss LSM 510 Meta confocal microscope.

For the assessment of Sox9 and laminin expression in wild-type and AmhCre Sox9 floxflox gonads at E13.5, embryos were fixed overnight in 4% paraformaldehyde (PFA) at 4 °C, then washed three times in 1× PBS. The embryos were processed and embedded into paraffin, cut at 5 μm, and mounted onto slides. The slides were baked at 60 °C (30 min), deparaffinized using three washes of xylene, and hydrated using three washes of 100% ethanol, then distilled water and 1× PBS. Antigen retrieval was performed by microwaving slides (on high) in 10 mM sodium citrate (pH = 6.0) for 20 min. Sections were then blocked for 30 min with 5% normal donkey serum, and stained overnight at 4 °C with primary antibodies for anti-Sox9 rabbit polyclonal (1:400) and anti-Laminin rabbit polyclonal (1:100). Sections were washed three times in 1× PBS with 0.1% Tween20 (1× PBST) and incubated with the fluorescent-conjugated secondary antibodies, Donkey anti-rabbit AlexaFluor488 (Thermo Fisher, 1 μg/mL), for 1 h at room temp. Sections were washed three times in 1× PBST, then incubated in 0.1% Sudan Black in 70% EtOH for 5 min to quench background autofluorescence. Lastly, sections were washed three times in 1× PBST, counterstained using DAPI, then washed three times in 1× PBS and mounted using Dako Fluorescent Mounting Medium (Dako). Sections were imaged using fluorescence microscopy (Olympus Corp).

46,XY DSD cases with uninformative exome sequencing

C57BL/6J-Y poschiavinus mice as a model for 46,XY gonadal dysgenesis

In cases where no pathogenic variant was found by exome sequencing, we identified many VUS outside of the DSD clinical gene list. To investigate the relevance of these VUS in regard to patients’ phenotype, we utilized a powerful mouse model for studying undervirilization in human 46,XY individuals. In this model, the presence of a Y chromosome originating from a M. domesticus poschiavinus strain (Y ^POS ) on a C57BL/6J (B6) background (B6-Y ^POS ), an inbred laboratory strain that normally carries a M. musculus Y chromosome, results in disrupted testicular development and female genital phenotype [38].

This inherited phenomenon has been extensively studied in the B6-Y ^POS animals. The failure to develop testes stems from the inability of the Sry ^POS gene to initiate normal testicular development when B6 autosomal and/or X-linked factors are present. Virtually all B6-Y ^POS animals develop some ovarian tissue; half develop exclusively ovarian tissue, classified as completely sex-reversed; and the remainder develop both ovarian and testicular tissue, classified as partially sex-reversed (gonad morphology shown in Fig. 1a). The B6-Y ^POS mice represent a good model for studying 46,XY DSD with gonadal dysgenesis because of the overlap of the major phenotypic features that are present in both humans and mice such as normal physical appearance without clinical findings including major organs other than the reproductive system, normal karyotype, external genitalia that range from ambiguous to typical female-like, internal genitalia ranging from absent Müllerian structures to presence of a uterus, and abnormal gonadal development characterized as dysgenetic testes, streak, or ovotestes.

In the embryonic mouse gonad, Sry is normally expressed in a dynamic wave (central to distal) between E10.5 and E12.5 in the XY genital ridge, with peak Sry expression occurring in normal XY ^B6 genital ridges at ~E11.5, i.e., at the 16–18 tail somite stage of development, which is followed by the upregulation of Sox9 [39, 40]. In contrast, expression of the Sry ^POS gene peaks 10 to 14 h later in the genital ridges of B6-Y ^POS fetuses [41]. We hypothesized that abnormal gonadal expression of specific genes in B6-Y ^POS males, after the surge of Sry during gonadal development, would correlate with the genes in which VUS were identified in 46,XY DSD patients by exome sequencing.

Gene expression differences between B6-Y ^B6 and of B6-Y ^POS males

Since all of the 46,XY DSD patients in the cohort carried a functional SRY gene, it was important to perform gene expression analysis in the animal model after the peak of Sry expression for optimal comparability. To achieve this, gonadal tissue from WT B6-Y ^B6 and undervirilized B6-Y ^POS males at embryonic day E11.5, specifically at 21 tail somites (a time point when the surge of Sry gene was complete in both B6-Y ^B6 and B6-Y ^POS males), was collected to perform RNA sequencing for assessment of differential gene expression.

Using this method, we identified 515 genes that were differentially expressed between B6-Y ^B6 and B6-Y ^POS males with a fold change greater than 1.5. Out of these 515 genes, 308 were underexpressed and 207 were overexpressed in B6-Y ^POS males (Fig. 1b; Additional file 2: Table S2). To validate the integrity of the method and make sure that the correct tissue was dissected at the correct embryonic stages, we first looked at the expression levels of two important genes involved in testicular development, Sry and Sox9. High Sry and low Sox9 expression levels in B6-Y ^POS males indicated the correct timing of embryonic development (Fig. 1c), expression of which coincided with a previous publication [41]. Second, we verified which genes in our DSD gene list used for exome variant filtering were present in the B6-Y ^B6 /B6-Y ^POS differentially expressed gene list. The comparison of the two lists revealed that 21 genes were in common: 15 underexpressed and 6 overexpressed (Additional file 2: Table S2). In our previous cohort [22], out of these 21 genes, 3 (HSD17B3 (hydroxysteroid 17-beta dehydrogenase 3), STAR (steroidogenic acute regulatory protein), FGFR2) contained a pathogenic variant identified by exome sequencing, explaining a total of 5 cases, and 2 (DHH (desert hedgehog), MAMLD1 (mastermind-like domain-containing 1)) contained a variant that was reported to the clinician to orient further endocrine or imaging testing toward a definitive diagnosis. Cumulatively, these findings indicate that the gene expression analyses were carried out in a correct tissue type, at the correct developmental time point, and that the differentially expressed genes between B6-Y ^B6 and B6-Y ^POS males may potentially be important in sex development.

Filtering of VUS in 46,XY DSD cases using the B6-Y ^POS gene list

On average, exome sequencing identifies ~ 21,000 variants per single case [23]. Since DSDs are rare conditions, all variants identified in exome with a minor allele frequency (MAF) of more than 1% in the population were excluded. The variants remaining after the MAF cutoff were classified as variants of unknown significance. The number of genes with a VUS in each case ranged from 30 to 1100 with an average of approximately 730 genes per case. The gene list generated via expression studies in B6-Y ^B6 /B6-Y ^POS mice, consisting of 515 genes, was used to filter the list of VUS-containing human genes identified by exome sequencing. The comparison of two lists identified 305 (189 underexpressed and 116 overexpressed in B6-Y ^POS ) genes that were both differentially expressed in B6-Y ^POS males and contained a VUS in the 46,XY DSD cohort with an uninformative exome (Fig. 1b).

All these genes are known to be expressed in the developing gonad at the time of sex determination (e.g., the method used to identify these genes intrinsically already ensures that all those genes are expressed in the relevant tissue (gonad) at the relevant developmental time). In order to increase the probability of identifying relevant candidate genes involved in 46,XY DSD pathogenesis, we further queried if the differentially expressed genes from the mouse model (all 515 genes) were involved in any known biological processes. Gene Ontology Consortium (GOC) [42] enrichment analysis confirmed that genes underexpressed in B6-Y ^POS males were indeed enriched in biological processes known to control multicellular organism and anatomical structure development, including male reproductive development (Additional file 3: Table S3). Understanding the relevance of the genes that were overexpressed in B6-Y ^POS males was less straightforward. These genes were enriched in only two biological processes: response to extracellular stimulus and epithelial cell differentiation. Both of these categories had a high P value indicating that many genes in the overexpressed category are not associated with any known biological processes at this time. In addition, all of the pathogenic variants identified in our previous 46,XY DSD cohort [22] were in the underexpressed category of genes, indicating that they need to be expressed at higher levels in the developing gonad for proper testicular formation.

Based on these findings, we focused on variants identified in genes underexpressed in B6-Y ^POS males whose higher expression in WT males correlated with normal male sex development. To choose a fold change cutoff, we looked at fold change differences in expression between B6-Y ^B6 and B6-Y ^POS males for genes present in our clinical primary gene list. We found that the majority of the genes have an expression that is 2-fold higher in WT males compared in B6-Y ^POS males (Fig. 1d). To make the analysis more stringent and improve the confidence of identifying true candidate genes involved in male sex development, we therefore increased the fold change cutoff in expression between the B6-Y ^B6 and B6-Y ^POS males from 1.5 to 2. This change decreased the number of underexpressed genes from 189 to 53. Additional filtering was performed based on variant frequency (variants with MAF close to or below 0.1%), amino acid conservation (variants in highly conserved residues across multiple species were given preference), number of variants contained in a gene across the cohort, in silico predictions for pathogenicity (preference was given to the variants predicted to be deleterious or damaging), availability of literature (some weight was given to genes with known functions), and gonadal cell-specific expressivity using GenitoUrinary Developmental Molecular Anatomy Project (GUDMAP) data [43] (preference was given to genes expressed in male-typical cells).

Using the abovementioned filtering criteria, we identified 15 novel candidate sex developmental genes, variants in which may be involved in 46,XY DSD pathogenesis. The list of VUS identified in the 46,XY cohort is shown in Table 2. The relative expressions of these genes in B6-Y ^B6 , B6-Y ^POS , and WT females are shown in Fig. 2a. The expression changes of all 15 genes were confirmed using quantitative real-time PCR (Fig. 2b).

Gene	DSD case ID	Zygosity	HGVSc	HGVSp	MAF gnomAD (%)
TOX2	RDSD021	Het	c.319G>A	p.Gly107Ser	0
	CDSD036	Cmpd Het	c.448A>G	p.Met150Val	0
	CDSD036	Cmpd Het	c.1201C>G	p.Pro401Ala	0
	CDSD036	Cmpd Het	c.1122_1124dupGCC	p.Pro376dup	0
DUSP15	RDSD020	Het	c.563G>C	p.Arg188Pro	0.002
NKD2	RDSD003	Het	c.1151G>A	p.Arg384Gln	0.001
CNGA1	CDSD030	Het	c.1478G>A	p.Arg493Gln	0.09
	RDSD022	Het	c.398G>T	p.Gly133Val	0.03
PTK2B	RDSD011	Het	c.1799G>A	p.Arg600Gln	0.0008
ESPN	RDSD044	Het	c.2230G>A	p.Asp744Asn	0.02
SMOC2	CDSD030	Het	c.1276G>A	p.Val426Met	0.3
ADAMTS16	RDSD013	Het	c.2200G>A	p.Val734Ile	0.8
	RDSD002	Het	c.298C>T	p.Arg100Trp	0.1
	RDSD022	Het	c.1405T>G	p.Phe469Val	0.02
FBLN2	CDSD030	Het	c.1486G>A	p.Ala496Thr	0.033
	CDSD029	Het	c.3605C>G	p.Ala1202Gly	0.004
NIPAL1	RDSD003	Het	c.1207A>G	p.Thr403Ala	0.1
	CDSD031	Het	c.31G>A	p.Glu11Lys	0
CYP26B1	CDSD032	Het	c.805C>G	p.Leu269Val	0.008
SPRY4	CDSD039	Het	c.446C>G	p.Pro149Arg	0.0004
MYBL1	RDSD004	Het	c.754T>A	p.Phe252Ile	0.05
	CDSD029	Het	c.1832G>C	p.Ser611Thr	0.0008
	RDSD049	Het	c.936T>A	p.Asn312Lys	0.03
ETV4	RDSD006	Het	c.523C>A	p.His175Asn	0.1
LGR5	RDSD007	Het	c.1834G>A	p.Val612Met	0.004
	RDSD020	Het	c.2341C>G	p.Pro781Ala	0.8
	RDSD048	Het	c.2537C>A	p.Thr846Asn	0

Het heterozygous, Cmpd het compound heterozygous, HGVSc Human Genome Variation Society coding sequence location, HGVSp Human Genome Variation Society protein sequence location, MAF minor allele frequency, gnomAD genome Aggregation Database

Expression of the novel candidate genes is Sox9-dependent

The time point chosen for our gene expression analysis was such that the Sry gene expression was similar between B6-Y ^B6 and B6-Y ^POS males. At that time, the downstream target of Sry, Sox9 was significantly decreased in B6-Y ^POS males (Fig. 1c), as previously described [41]. In order to identify if Sox9 had any effect on expression of the candidate genes, we used the Amh-Cre Sox9flox/flox mouse model where Sox9 expression is suppressed in Sertoli cells [44]. By E13.5, Sox9 protein is completely absent (Fig. 3a), and these mice show postnatal fertility defects [44]. Earlier Sox9 knockout models result in complete sex reversal (XY with ovaries) [45, 46] or embryonic lethality [47]; neither situation sheds light on Sox9 target genes during sex determination. The Amh-Cre Sox9flox/flox mouse model allows the examination of Sox9 loss in an intact Sertoli cell environment.

Performing gene expression analysis via RNA sequencing in mice with suppressed Sox9 expression showed that 13 of the novel candidate genes for 46,XY DSD underexpressed in B6-Y ^POS males, were also underexpressed in Sox9 knockout male gonads with 7 being significantly different (Fig. 3b). This finding suggests that Sox9 may be upstream of some of the novel candidate genes for 46,XY DSD. In addition, the profiles of gonadal gene expressions from GUDMAP reveal that in almost all cases, the patterns of gene expression are similar to bona fide target genes of Sox9 such as Amh (anti-Müllerian hormone) and Ptgds (prostaglandin D2 synthase) [42, 48, 49]. The target gene expression is higher in the male supporting cells (Sertoli) than in the female supporting cells (granulosa) (Fig. 4). The rest of the genes may be regulated by other transcription factors such as Nr5a1, which is a known regulator of Cyp26b1 (cytochrome P450 family 26 subfamily B member 1) [50]. Collectively, our results show that variants in the candidate genes such as the ones we have identified in the 46,XY DSD cases (Table 2) may be responsible for the patient’s phenotype.