Comprehensive molecular analysis identifies eight novel variants in XY females with disorders of sex development

Abstract

Disorders of sex development (DSD) are a group of clinical conditions with variable presentation and genetic background. Females with or without development of secondary sexual characters and presenting with primary amenorrhea (PA) and a 46,XY karyotype are one of the classified groups in DSD. In this study, we aimed to determine the genetic mutations in 25 females with PA and a 46,XY karyotype to show correlations with their phenotypes. Routine Sanger sequencing with candidate genes like SRY, AR, SRD5A2, and SF1, which are mainly responsible for 46,XY DSD in adolescent females, was performed. In a cohort of 25 patients of PA with 46,XY DSD, where routine Sanger sequencing failed to detect the mutations, next-generation sequencing of a targeted gene panel with 81 genes was used for the molecular diagnosis. The targeted sequencing identified a total of 21 mutations including 8 novel variants in 20 out of 25 patients with DSD. The most frequently identified mutations in our series were in AR (36%), followed by SRD5A2 (20%), SF1 (12%), DHX37 (4%), HSD17B3 (4%), and DMRT2 (4%). We could not find any mutation in the DSD-related genes in five (20%) patients due to complex molecular mechanisms in 46,XY DSD, highlighting the possibility of new DSD genes which are yet to be discovered in these disorders. In conclusion, genetic testing, including cytogenetics and molecular genetics, is important for the diagnosis and management of 46,XY DSD cases.

Introduction

The disorders of sex development (DSD) are congenital conditions in which development of chromosomal, gonadal, or anatomical sex is atypical (Strauss et al., 2017). Adolescent females with DSD and presenting with primary amenorrhea (PA) have variable phenotypic presentations depending on the defects of gonadal development (testis/ovary) (e.g. complete gonadal dysgenesis known as Sweyer syndrome) or conditions that lead to defective hormone biosynthesis or action of androgens (e.g. androgen insensitivity syndrome, AIS) (Buonocore et al., 2019). Among DSD, the incidence of phenotypic females with a 46,XY karyotype and presenting with PA is about 14%, as estimated by the practice committee of American Society of Reproductive Medicine (Practice Committee of the American Society for Reproductive Medicine, 2008). These conditions are estimated to occur with a frequency of 1:10 000 individuals in the population (Wang et al., 2018).

Development of sex is mainly dependent on genetic sex, which in turn decides the gonadal sex, and as genital development goes by way of the genetic and gonadal sex, the phenotypic sex of the individual is determined. Several genes are involved in the complex process of sex determination in humans and malfunctioning of any one of them leads to DSD (Paliwal et al., 2011). Development of gonads is bipotential during the first 7–8 weeks of mammalian embryogenesis, and subsequently differentiation occurs to give rise to testes or ovaries in the individuals with a 46,XY or 46,XX karyotype, respectively (Paliwal et al., 2011). The individuals with normal 46,XX (female) karyotype and normal ovaries should develop female external and internal genitalia. But, there are individuals who have 46,XY (male) karyotype with testes or dysgenetic gonads and labeled as female or ambiguous genitalia at the time of birth, and subsequently grow up as females. Such females present with PA with or without development of normal secondary sexual characters at the time of puberty. This can be a social stigma for the individual as well as the family members and may cause a psychological burden for them. Despite having the Y chromosome, such females may not develop normal testes and normal male internal and external genitalia. The reason for non-development of testes is a defect in the genes that determine the sex of these patients, due to which the bipotential gonads fail to develop into either testes or ovaries. The sex-determining genes are expressed in the developing genital ridges and the differentiation of gonads into one sex or the other is a complex process due to the interaction of the various gene products responsible for gonadal and genital differentiation (MacLaughlin and Donahoe, 2004). If the phenotypic females develop normal testes, the location of such testicular tissue is at an abnormal region and there may be defective androgen production or failure of conversion of testosterone into its active form, dihydrotestosterone (DHT), by the enzyme steroid 5 alpha-reductase (SRD5A2 gene defects) or, more often, there may be defective action of androgen due to defective androgen receptors (AR gene defects). This leads to a problem with masculinization of the external genitalia giving rise to the abnormal phenotype appearing as female (Practice Committee of the American Society for Reproductive Medicine, 2008). Thus, the majority of the females with an XY karyotype and PA are diagnosed as having AIS (Minto et al., 2003).

However, revelation of molecular mechanisms behind the other potential causes of PA could be possible due to recent advancements like next-generation sequencing (NGS) (Maimoun et al., 2011). A genetic diagnostic rate of 43% has been reported by using a targeted DSD gene panel in a large patient cohort of 278 patients with 46,XY DSD (Eggers et al., 2016). In view of the importance of the genes responsible for DSD, this study was undertaken with an aim to identify the underlying molecular changes by using NGS of a targeted gene panel to investigate a cohort of 25 patients with PA and 46,XY DSD.

Materials and methods

Ethical statement

This study was approved by the Institutional Ethics Committee for research on human patients of the ICMR-National Institute of Immunohaematology, Parel, Mumbai, India.

Clinical assessment of the patients

We have studied 25 patients of PA with 46,XY DSD and no menarche by the age of 13 or later, regardless of the presence of normal growth and development, with or without the appearance of the secondary sexual characters. The patients’ details like age, height, clinical features, breast development, pubic hairs, axillary hairs, hormone profile, and ultrasound details were recorded in the case record sheet. Detailed family histories were collected, and pedigree charts were drawn. Pedigree charts of the patients with novel variants and patients with previously reported mutations are shown in Supplementary Figs S1A, B, C, D, E, F, G, and H and S2A, B, C, D, E, F, G, H, I, J, K, and L, respectively. The study protocols were approved by the Institutional Ethics committee for human patients. Informed consent was obtained from all adult participants and in case of adolescents, the consent was obtained from their parents prior to peripheral blood collection.

Conventional cytogenetics and FISH

The whole peripheral blood collected in heparin vacutainers was cultured (37°C for 72 h) in RPMI 1640 medium (9 ml) (Gibco, New York, USA), supplemented with bovine serum (1 ml) (Gibco), l-glutamine (0.1 ml), phytohemagglutinin (Gibco), and antibiotics (penicillin and streptomycin; Gibco). The cultures were treated with a hypotonic solution (0.075 M KCl), fixed in methanol:acetic acid (3:1), and dropped on pre-chilled slides. The chromosomal preparation was subjected to GTG banding and karyotyped according to International System of Human Cytogenetic Nomenclature (ISCN 2016) (Moorhead et al., 1960; Seabright, 1971).

The chromosomal preparations were passed through a series of alcohol (70%, 80%, and 100%) for 5 min at room temperature and then air dried. The centromeric and locus-specific probes for the X and Y chromosomes (Vysis, Abbott, USA) were used for FISH, and hybridization was carried out using a standard protocol (Davies, 1993). The FISH probes were prepared as per the manufacturer’s instructions and applied. The co-denaturation of the prepared slides was performed at 73°C for 4 min. After co-denaturation, slides were covered with coverslips by applying rubber cement and kept for hybridization at 37°C for 24 h. The post-hybridization washing was performed using wash solution of different stringency (low and high) at 72°C. The slides were counterstained with DAPI (4′,6-diamidino-2-phenylindole, 7.5 μl/slide; Sigma-Aldrich, St. Louis, MO, USA), and analyzed under a fluorescent microscope (Nikon 90i). The images were captured with a CCD camera and analyzed with Applied Spectral Image System (ASI, Israel).

Molecular study

Genomic DNA was extracted from peripheral blood using a human blood DNA extraction kit (Qiagen GmbH, Hilden, Germany) as per the standard protocol. Sanger sequencing was carried out for SRY, AR, SRD5A2, and SF1(NR5A1) genes and primers were designed using Primer 3 software (Supplementary Table SI). The annealing temperatures for the amplification of the SRY, AR, SRD5A2, and SF1 genes are shown in Supplementary Table SI. Sequencing was carried out using an ABI-3130 genetic analyzer, the analysis was performed using Chromas Lite software (Technelysium Pty Ltd, South Brisbane, Australia), and nucleotide sequences were compared with the wild-type DNA sequences.

Next-generation sequencing

NGS was carried out for a DSD gene panel covering 81 DSD-related genes (Supplementary Table SII). We performed targeted exome sequencing at Med-Genome Labs Pvt. Ltd. (Bangalore, India) using the custom capture kit. Sequencing (Illumina platform, San Diego, CA, USA) with a mean coverage of 80–100× were performed. Using the BWA program, the sequences were then aligned to the human reference genome (GRCh38.p13) and analyzed using Picard and GATK version 3.6. The VEP program was used to annotate the gene variants. Furthermore, the clinically relevant mutations were annotated using datasets such as ClinVar, OMIM (update on 20 February 2020), GWAS, HGMD (v2019.4), and SwissVar. Furthermore, the common variants were filtered based on allele frequency in 1000 Genome Phase 3, gnomAD (v2.1), EVS, dbSNP (v151), 1000 Genome, and an internal Indian population database.

Bioinformatics analysis

Primer 3 software (online available) version 4.1.0 was used to design the primers (https://primer3.ut.ee/). SRY, AR, SRD5A2, and SF1 gene references were obtained from the Ensembl genome browser (https://asia.ensembl.org/index.html) and NCBI (https://www.ncbi.nlm.nih.gov/). Sequencing data were analyzed with a reference genome for screening the mutations using the Clustal Omega tool (https://www.ebi.ac.uk/Tools/msa/clustalo/) (Hinxton, Cambridgeshire, UK). The pathogenicity of the novel missense mutations was confirmed using online bioinformatics tools such as Polyphen-2 (http://genetics.bwh.harvard.edu/pph2/), SIFT (http://sift.jcvi.org/), and MutationTaster (www.mutationtaster.org/). The variants were further classified based on ACMG guidelines (Richards et al., 2015).

Protein modeling

The given protein sequences of AR, SRD5A2, DMRT2, and SF1/NR5A1 were considered for model building using different software based on their sequence identity and query coverage. To begin with, the protein sequence of the wild-type SRD5A2 was modeled using a homology modeling-based method on the SWISSMODEL online server (Guex et al., 2009). To understand the impact of the mutation on the protein structure, the mutant protein was also modeled using the same server. DMRT2, AR, and SF1/NR5A1 were modeled using the Robetta online server, which is a deep learning-based method for model building (Kim et al., 2004). All of the generated models were subjected to energy minimization using the GROMACS standalone server (Berendsen et al., 1995), wherein topology files were generated using OPLS AA/L force field. All of the 3D models were subjected to 50 000 steps during energy minimization. Next, these optimized models were subjected to structure validation using ProSA-Web (Wiederstein and Sippl, 2007) and PROCHECK-Ramachandran plot (Laskowski et al., 1993) online servers, for the local energy and phi-psi angle analysis, respectively.

Protein–ligand docking

The modeled SRD5A2 protein was examined for docking with NADPH using Hex standalone software (Ghoorah et al., 2013). The 3D ligand of the NADPH was extracted from PDB ID: 7BW1 (Xiao et al., 2020). Regarding docking, the wild-type SRD5A2 protein was loaded in to HEX software. Within the software, ‘correlation type’ was set as shape + Electro + DARS. Furthermore, ‘Post-Processing’ was DARS minimization. The grid dimension was maintained at 0.6 Å. Finally, post-docking generated 10 poses that were subjected to CHIMERA (Pettersen et al., 2004) for structure visualization and LIGPLOT (Laskowski and Swindells, 2011) for amino acid interaction analysis. Similar steps were followed for the docking of the mutant protein with the NADPH ligand.

Results

Clinical evaluation

The study included 25 unrelated patients with DSD who mainly presented as PA. The common clinical features associated with these cases were absence of secondary sexual characters, ambiguous genitalia, absent Müllerian structures, and gonadal abnormality in the form of either absent ovaries or streak gonads or normal or undescended testis in phenotypic females (Table I). The pedigree analysis of all the patients showed that there was no history of parental consanguinity except for Patients P13 (Supplementary Fig. S2H) and P20 (Supplementary Fig. S1H). None of the parents or siblings had similar complaints to that of the 46,XY DSD patients except for Patients P4 (Supplementary Fig. S1B) and P5 (Supplementary Fig. S2C), where siblings were also affected. The age ranged from 13 to 31 years with a mean age of 19.2 years. The clinical examination revealed normal female external genitalia in 56% of cases and ambiguous genitalia (pseudo-vaginal perineal hypospadias or clitoromegaly) in 44% of cases (Fig. 1). Clinical evaluation of breast development in XY DSD cases with PA showed Tanner stage 1 (48%), stage 2 (20%) and stage 3 (32%), with no patients showing stage 4 or 5. The uterus was found to be absent in 76% of cases. Gonads were found to be testis in 44% of cases, streak gonads in 24% of cases, and absent in 32% of cases (Fig. 1). The levels of FSH and testosterone were recorded and are shown in Table I.

Cytogenetic analysis

Chromosomal analysis revealed a 46,XY karyotype in each of the 25 patients. A FISH study showed a XY karyotype in 22 cases and mosaic patterns for 2 cases with X/XY [46,XY (90%), 45,X (10%)] and 1 case with XX/XY [46,XX (60%), 46,XY (40%)].

Molecular study

The targeted sequencing identified a total of 21 mutations in 20 out of 25 patients with DSD. The overall incidence was found to be 84%. The most frequently identified mutations in our series were in the AR gene (nine patients, 36%), followed by SRD5A2 (five patients, 20%), SF1 (three patients, 12%), DHX37 (one patient, 4%), HSD17B3 (one patient, 4%), and DMRT2 genes (one patient, 4%) (Fig. 2). In our study, we found that 8 of the 21 mutations (38.09%) were novel variants, of which three were in the AR gene, two were in the SRD5A2 gene, and one each was in the SF1, DMRT2 and HSD17B3 genes.

In the present study, we found three mutations, c.175C>T (p.Gln59Ter), c.1443C>G (p.Tyr481Ter), and c.1567G>T, (p.Glu523Ter) in Exon 1 of the AR gene, two missense mutations, c.1742A>C (p.Lys581Thr) and c.1762G>C (p.Ala588Pro) in Exon 2 of the AR gene, and three mutations, c.2226G>A (p.Trp742Ter), c.2255G>A (p.Trp752Ter), and c.2301del (p.Asp768IlefsTer21) in Exon 5 of the AR gene, while one patient had a missense mutation c.2323C>T (p.Arg775Cys) in Exon 6 of the AR gene (Table II).

The SRD5A2 gene was found to be the second most frequently mutated gene in our study associated with DSD. A total of six mutations in SRD5A2 were identified in five patients. Out of these, three were in Exon 1, two were in Exon 4, and one was in Exon 5. The mutation c.169G>T (p.Glu57Ter) was identified in two cases. The mutations c.81_94del, c.589G>A (p.Glu197Lys), c.691C>T (p.His231Tyr) and c.737G>A (p.Arg246Gln) were identified in one case each (Table II).

Three mutations were identified in the SF1 (NR5A1) gene; c.19G>T (p.Glu7Ter), c.97T>C (p.Cys33Arg) in Exon 2 and c.250C>T (p.Arg84Cys) in Exon 4.

One mutation was found in each of the DHX37, DMRT2, and HSD17B3 genes; c.1877C>T (p.Ser626Leu), c.71A>C (p.Glu24Ala), and c.385 + 5G>A, respectively (Table II).

The electropherograms of Sanger sequencing of all the novel and reported mutations are shown in Fig. 3 and Supplementary Figs S3 and S4, respectively.

Protein modeling of the novel variants

DMRT2 p.Glu24Ala

The 3D structure of the DMRT2 protein was generated using the Robetta server and considered for hydrogen bond analysis using CHIMERA standalone software. In particular, Glu24 of the wild-type DMRT2 protein shows a salt bridge formation with Arg31 (Fig. 4A). However, after the p.Glu24Ala substitution, because the methyl side chain of alanine is non-reactive, there is no salt bridge formation with Arg31 (Fig. 4B). Thus, the mutation results in the loss of salt bridge formation, which could certainly have a local impact on the flexibility of the loop region of the protein. A literature review confirmed that any substitution of charged residue with alanine has an adverse effect on protein stability (Pakula et al., 1986; Cunningham and Wells, 1989; Yang et al., 2009).

AR p.Glu523Ter, AR p.Trp742Ter

The 3D model of the protein sequence of AR showed potential termination residues at positions like 523, 742, and 752. Position 523 harbors glutamic acid whereas Positions 742 and 752 have tryptophans. Mutations at Positions 523, 742, and 752 introducing stop codons could terminate transcription of the ligand-binding domain (LBD). Thus, truncation removes the ligand-binding helical secondary structure, which is critically involved in receptor dimerization (Xiao et al., 2020) (Fig. 4C).

AR p.Lys581Thr

Another 3D model of the protein sequence of AR in wild-type form showed two interactions at Lys581 where there was a single bond formation with Cys 577 and Lys585 (Fig. 4D). After the mutation substituting threonine, there are three interactions observed wherein Thr581 is showing two interactions with Cys577 and a single interaction with Lys585 (Fig. 4E). This mutation could also be associated with the loss of an androgen binding site.

SF1 p.Cys33Arg

The generated wild-type model of SF1/NR5A1 has cysteine at Position 33. As a result, it forms a disulphide bridge with Cys30 and an additional interaction with Phe37. After the substitution of Cys with Arg, there is a loss of disulphide bridge formation (Fig. 4F and G). However, the hydrogen bond formation is still observed for Cys33 with Phe37, Cys30, and Val15.

SRD5A2 p.Ala28LeufsTer103

The given protein sequence of SRD5A2 was submitted to the SWISSMODEL server. Based on the template search, PDB ID: 7BW1 was listed as a potential 3D structure with a query coverage of 96% and sequence identity of 99.6%. This template was used for the modeling of the wild type and the mutant SRD5A2 protein. The SRD5A2 gene has a reported mutation wherein Ala28 is substituted by Leu due to the deletion of bases between 81 and 94. As a result, there is a frameshift from amino acid Position 29 and formation of a stop codon at Position 103, which results in a premature protein truncation. The frame shift changes the transmembrane and beta-strand region formation after the 29th position (Fig. 4H).

Furthermore, the structure validation report of all the modeled structures from PROCHECK and ProSA-Web report confirmed energetically stable structures.

SRD5A2 p.His231Tyr

Protein–ligand docking showed that NADP was docked within the active site pocket based on the template of 7BW1. When wild-type SRD5A2 was docked with NADP, the binding energy was −730.10 kcal/mol. The mutant with NADP showed a binding energy of −717.16 kcal/mol. Thus, a lower binding energy was observed after the substitution of histidine with tyrosine. Moreover, the histidine residue at Position 231 showcased hydrogen bond formation with NADP. However, after mutation, this hydrogen bonding ceases as per their hydrogen bond report (Fig. 4I and J). Based on the LIGPLOT report, there is an absence of hydrogen bonding between Tyr and NADP. Instead, there is an external bond formation observed with a light blue color (Fig. 5A and B).

HSDB17 intron 4′5 splice site

The HSDB17 mutation was not investigated by protein modeling as it created a change in one of the mRNA splice sites.

Discussion

DSD are rare congenital disorders, where there is discordance between chromosomal, gonadal, and the phenotypic sex of an individual. Although DSD are a major disorder of pediatric concern, many cases remain undiagnosed until puberty and the majority present as PA at the time of puberty (Mendonca et al., 2009). This study included a total of 25 patients, presumed to be females, with 46,XY DSD, who presented with PA with or without the development of secondary sexual characters and varying degrees of virilization. To the best of our knowledge, this is the first Indian study on such a large number of cases. We have systematically analyzed the genotypes of the cases and related them with their phenotypes, and we note the extreme variability in the clinical spectrum and genetic background of these cases. In the present study, we found that the phenotype of the DSD patients with PA is mainly female external genitalia in 56% of patients, while ambiguous genitalia in the form of pseudo-vaginal perineal hypospadias or clitoromegaly was also observed in 44% of cases, making the exact clinical diagnosis of these conditions difficult. Clinical evaluation of internal genitalia in these cases showed that majority of these cases (76%) were without Müllerian structures, and many cases (44%) had testes as gonads. The most common cause of the PA with DSD in this cohort was found to be AIS (36%) followed by steroid 5 alpha-reductase deficiency (25%), both of which result in defective androgen activity, and this was similarly found in a previous report (Costagliola et al., 2021). The cytogenetic investigations play an important role in establishing the diagnosis of DSD and should be advised immediately to those patients who are suspected to be having clinical features of 46,XY DSD. This is important, as post-pubertal complete AIS patients are prone to tumor development in undescended abdominal testes (Kathrins and Kolon, 2016). As per the literature, the estimated risk of malignancy is 15–50% in gonadal dysgenesis and <1–15% in disorders of androgen synthesis or action (Wisniewski et al., 2019). The 46,XY karyotype in females with PA and DSD have been reported with various frequencies (2–31%) in different populations (Ghosh et al., 2018). In the present study, in all the 25 DSD patients, the karyotype was found to be 46,XY with no additional chromosomal abnormalities. The chromosomal abnormalities using a combination of conventional G-banded cytogenetics and FISH revealed that 88% of the cases with PA (22/25) had a pure XY karyotype and the remaining 12% (3/25) had a mosaic karyotype. The frequencies of chromosomal abnormalities identified in our study were similar to the frequencies reported in the literature (García-Acero et al., 2020). Hence, the combination of cytogenetics tools is important in the accurate diagnosis of the sex abnormalities.

Various studies have been carried out in the past using targeted NGS, where mutations have been reported with different frequencies, commonly in AR, SRD5A2, SF1, DHX37, HSD17B3, SRY, DMRT1, DMRT2, and other DSD-related genes with low frequency (Wang et al., 2018; Costagliola et al., 2021; Yu et al., 2021). Our study is the single largest study, where we have evaluated molecular genetic characterization of PA cases with DSD, providing 80% of these cases with an exact molecular diagnosis with a targeted NGS approach among the Indian population. Reaching a specific diagnosis is important, not only to know the pattern of inheritance and chances of other family members being affected, but also for identifying the associated features and to know the risk of tumors in the long term (Achermann et al., 2015). The overall incidence of mutations in the AR, SRD5A2, and SF1 genes in our study was 72% (18/25) and was found to be similar as that of the study reported by Yu et al. (78.3%, 47/60). We have identified 8 novel variants out of the 21 variants reported (38.1%) in our study, which is similar to a frequency (41.7%, 25/60) previously reported in the literature (Yu et al., 2021). Our extensive literature search revealed that AR gene mutations ranging from a single nucleotide variation to complete gene deletion including intronic mutations (Brinkmann et al., 1996; Boehmer et al., 2001) are responsible for variable phenotypes in 46,XY DSD patients, leading to either complete or partial AIS (Nagaraja et al., 2019). The androgen receptor gene is present at the chromosomal location Xq11-12 and contains a total of eight exons. The AR gene encodes a protein of 919 amino acids, which is a member of the nuclear receptor super family and comprises an N-terminal domain (NTD) (Exon 1), a DNA-binding domain (DBD) (Exons 2 and 3), a hinge region (Exons 3 and 4), and an LBD (Exons 4–8). A defective AR protein can cause androgen insensitivity, known as AIS (OMIM # 300068) (Audí et al., 2010). We have identified nine cases (36%) with AR gene mutations, which confirms the diagnosis as AIS in these patients. So far, 32 different mutations in the AR gene have been reported in the Indian population (Nagaraja et al., 2019). Our single-centre study will increase these data with nine mutations in AIS patients. Out of the 32 mutations, most were found in the LBD region (14 mutations) and the least were found in the Hinge region (two mutations). Similarly, we also found most (4/9) of the mutations to be in the LBD region. In a study by Yu et al. (2021), the frequency of AR gene mutations was found to be 35% (21/60). We also identified the AR gene as the most common mutated (42.8%, 9/21) gene in our Indian population. In the present study, three novel variants in the AR gene were identified. The hemizygous nonsense variant c.1567G>T (p.Glu523Ter), found in Exon 1 of the AR gene, results in a stop codon and causes the premature truncation of the protein at Codon 523. Exon 1 of the AR gene encodes the NTD, and it mediates the transactivation function, the disruption of which compromises the biological function of the AR gene (Saranya et al., 2016). In silico prediction of the variant was found to be damaging by MutationTaster2 and it was found to be a pathogenic variant as per the ACMG guidelines. The second novel variant observed in our study was the missense hemizygous mutation c.1742A>C (p.Lys581Thr) leading to the replacement of the amino acid lysine with threonine at Position 581 in Exon 2 of the AR gene. Exons 2 and 3 of the AR gene encode the DBD (residues 539–628), which contain two peculiar zinc finger modules necessary to bind the DNA and stabilize the protein (Gulía et al., 2018). Many studies have reported a relatively high number of mutations that can damage the structure and function of this domain, causing AIS. This mutation was found to be likely pathogenic and in silico prediction of the variant was probably damaging by polyphen-2 and damaging by SIFT. Another novel variant observed was the nonsense mutation c.2226G>A (p.Trp742Ter) in Exon 5 replacing tryptophan at the amino acid Position 742 by a stop codon, leading to the premature truncation of protein. This mutation was found to be pathogenic and damaging by MutationTaster2. Exon 5 of the AR gene encodes the LBD of the protein. Receptor dimerization of the AR protein can occur due to strong interactions between the LBDs and weak interactions between DBDs (Brinkmann et al., 1996).

The second gene associated with 46,XY DSD in our study was the steroid 5 alpha-reductase 2 gene (SRD5A2) (cytogenetic location: 2p23.1). This gene consists of five exons and is responsible for conversion of testosterone to its active form DHT. Mutations in the SRD5A2 gene lead to pseudo-vaginal perineoscrotal hypospadias (OMIM # 264600) (Andersson et al., 1991), a rare autosomal-recessive (both the parents are carriers) disorder. We found that 6 out of 21 (28.6%) mutations were in the SRD5A2 gene and these were the second most frequent in our cohort. In a review by Nagaraja et al. (2019), a total of 19 mutations were in the SRD5A2 gene from various studies among the Indian population and they were found to be the second most frequently reported among the genes causing 46,XY DSD. Yu et al. (2021) also reported SRD5A2 gene mutations to be the second most common mutation (21.7%, 13/60) in an Asian population. In our study, we found three mutations in Exon 1, two in Exon 4, and one in Exon 5. Maimoun et al. (2011) reported a predominance of mutations in Exon 1 (35.8%), Exon 4 (21.7%), and Exon 3 (11.3%), while Exon 5 (9.4%) was found to be in relatively preserved regions. Various Asian studies have revealed that the majority of mutations were confined to Exons 1 and 5 of the SRD5A2 gene (Thigpen et al., 1992; Nie et al., 2011; Shabir et al., 2015), and this was similarly observed in Indian studies, where most of the mutations were observed in Exons 1 and 5, indicating that these exons are the hot spot regions of the gene (Eunice et al., 2008; Sahu et al., 2009; Nagaraja et al., 2010). Furthermore, we have identified two novel variants in the SRD5A2 gene. The novel homozygous mutation found in Exon 1 of the SRD5A2 gene was c.81_94del, with substitution of the amino acid alanine at Position 28, with leucine causing a frameshift from amino acid Position 29 and formation of a stop codon at Position 103, which resulted in the premature truncation of the protein. A mutated enzyme of 103 residues is therefore formed instead of a protein of 245 amino acids. This variant was found to be pathogenic as per ACMG guidelines. Point mutations leading to the formation of premature stop codons have been described in previous studies in patients with 5 alpha-reductase Type 2 deficiency (Hiort et al., 2002; Vilchis et al., 2008), but there are very few cases reported with mutations in the coding region of the SRD5A2 gene leading to premature termination of the enzyme. The mutations in Exon 1 of SRD5A2 were observed in different ethnic and geographical backgrounds (Mazen et al., 2003; Maimoun et al., 2011), but there are very few reports described with the presence of heterozygous mutations (Thigpen et al., 1992).

Another novel mutation observed in our study is the compound heterozygous mutation c.691C>T in Exon 4 that replaces the amino acid histidine with tyrosine at amino acid Position 231. This mutation was found to be likely pathogenic as per the ACMG guidelines and was predicted to be probably damaging by Polyphen 2. The positively charged histidine at Position 231 shows hydrogen bond formation with NADP. However, after mutation, this hydrogen bonding ceases, as per the CHIMERA report. Based on the LIGPLOT report, similar to the CHIMERA report, the hydrogen bonding was absent between the tyrosine and NADP (Fig. 5A and B). Instead, there is an external bond which may disrupt the cofactor binding, which is quite evident from the mutation.

The other gene that is responsible for 46,XY DSD is SF1 (Steroidogenic Factor 1) (OMIM # 184757), located at 9q33, also known as the NR5A1 gene. It is an orphan nuclear receptor that consists of 461 amino acids and is structurally similar to the other nuclear receptors (Oba et al., 1996; Wong et al., 1996). It encodes a nuclear transcription factor regulating the expression of several genes that participate in sexual development. In humans, heterozygous SF1 mutations in XY individuals lead to adrenal and gonadal failure (Köhler and Achermann, 2010), cryptorchidism (Wada et al., 2005), micro-penis and infertility (Bashamboo et al., 2010). Fabbri-Scallet et al. (2018) reviewed 238 published cases of both 46,XX and XY DSD and observed SF1 gene mutations at a frequency of 12% (25/205) in 46,XY DSD cases, with a predominance of mutations in the DBD (35%) and LBD (42.3%) domains. The frequency of SF1 gene mutations leading to DSD in our study was found to be 14.3% (3/21), similar to the study by Fabbri-Scallet et al. We found three heterozygous mutations in the DBD encoding region of SF1: two in Exon 2 and one in Exon 4. Nagaraja et al. (2019) reviewed four cases of SF1 gene mutations in gonadal dysgenesis among the Indian population, out of which three were reported by in male patients and one was reported by Paliwal et al. (2011) in a female patient. We identified one novel variant c.97T>C (p.Cys33Arg), a likely pathogenic, missense mutation in the DBD region of the protein. In silico prediction of this variant is probably damaging by polyphen-2 and damaging by SIFT. Fabbri-Scarlet et al. (2018) reported a mutation on same position; however, the amino acid change reported was a serine instead of a cysteine.

The less frequently identified gene mutations observed in our study were in one patient each with DMRT2, DHX37, and HSD17B3 (4.76%, 1/21) mutations and diagnosis of 46,XY DSD. Previous studies also reported these mutations as less frequent mutations among 46,XY DSD females (Yu et al., 2021). We have identified a novel heterozygous missense variant in Exon 3 of the DMRT2 gene in a female with PA and it was not reported previously in 46,XY DSD cases. This variant resulted in the amino acid substitution of alanine for glutamic acid (E) at Codon 24 (p.Glu24Ala). In silico prediction of the variant was damaging by SIFT. Although this variant is classified as a variant of uncertain significance as per ACMG guidelines, parental testing is essential and may change the classification. Unfortunately, we could not perform the parental sequencing, as the parents were not available for the molecular study. The DMRT2 gene (encoding doublesex and mab-3-related transcription factor 2) has six exons and is believed to be involved in gonadal development. To date, only eight clinically significant DMRT2 gene mutations have been reported from ClinVar and Humasvar with links to dbSNP (https://www.genecards.org/).

The DEAH-box RNA helicase (DHX37), essential for ribosome biogenesis, is found specifically in association with 46,XY gonadal dysgenesis (OMIM # 273250) and 46,XY testicular regression syndrome (TRS). The pathogenic mutations in the DHX37 gene (OMIM # 617362) are responsible for phenotypes ranging from phenotypic females to males with bilateral or unilateral undescended testis. The exact prevalence of TRS is not known but it has been found to affect ∼1:2000 boys (Pirgon and Dündar, 2012). In our study, we found one patient with a heterozygous missense variation, c.1877C>T, in Exon 15 of the DHX37 gene that resulted in the amino acid substitution of leucine for serine at Codon 626 (p.Ser626Leu). A similar variation was reported by McElreavey et al. (2020) in a male baby with TRS.

Another novel intronic splice-site variant, c.385 + 5G>A, was found in the HSD17B3 gene in a single patient with DSD and PA in our study. Homozygous or compound heterozygous mutations in the HSD17B3 gene lead to 17-beta hydroxysteroid dehydrogenase III (17β-HSD3) deficiency (OMIM # 264300). In the largest DSD European database, patients with 17β-HSD3 deficiency represent ∼4% of total 46,XY DSD patients (Hughes, 2008). Mendonca et al. (2017), found 37 mutations in the HSD17B3 gene reported in the literature and they were found to be missense mutations, nonsense mutations, exonic deletions, duplications, intronic splice sites, and amplification mutations. This disorder is characterized by hypoplastic to normal internal genitalia (epididymis, vas deferens, seminal vesicles, and ejaculatory ducts) with female external genitalia and the absence of a prostate.

Genotype–phenotype correlation in patients with DSD and clinical relevance

It is difficult to establish a clear genotype-phenotype correlation in patients with 46,XY DSD as there are multiple genes involved in the etiopathogenesis of DSD (Maimoun et al., 2011). Those females who present with PA and karyotype 46,XY are phenotypically female since the abnormal gonadal tissue in these cases fails to produce Müllerian inhibiting factor and testosterone and, even if they produce testosterone, there is a defective androgen receptor for its action or it fails to be converted into the active form, DHT. The genotype–phenotype correlation was assessed in 25 unrelated patients with 46,XY DSD and presenting with PA with molecular defects in AR, SRD5A2, SF1, DHX37, DMRT2, and HSD17B3 genes (Tables I and II). All the patients from P1 to P9 with a molecular diagnosis of AR gene mutations were phenotypic females with 46,XY karyotypes (except Patient P3 who had a mosaic karyotype), with variable degrees of Tanner staging of breast, pubic and axillary hair development, absent uterus and cryptorchidism of one or both gonads. As all the patients with AR gene mutations, including the novel variants, have normal female phenotypes and presented at the time of puberty with PA, the diagnosis of partial androgen insensitivity syndrome (PAIS) is less likely. We were not able to establish the inheritance pattern of novel and reported mutations as the parents were not available or refused to participate in the study and this is the limitation of the present study. In Patients P1, P2, P4, and P9, testicular tissue was not located with ultrasound findings, needing repeat evaluation or confirmation by MRI. Elevated levels of testosterone were observed in three patients (P1, P2, and P8), which is typically seen in AIS. In the other six cases, although testosterone levels were not known to define a specific clinical diagnosis, molecular evaluation confirmed the diagnosis of AIS in these cases. The management of such cases should include creation of a functional vagina and removal of cryptorchid gonads as these may lead to development of malignancy. Gonadectomy can be delayed in those with AIS as gonadal tumors rarely occur before puberty in such patients and the limited pubertal development results from endogenous gonadal hormone production rather than exogenous hormone treatment (Speroff and Fritz, 2005). The cryptorchid testis has a relatively high incidence of neoplasia; therefore, a gonadectomy is suggested before exogenous hormone treatment. The prevalence of a germ cell tumor is 15% in PAIS and the risk increases after puberty, reaching up to 33% after the age of 50 (Araujo-Melendez et al., 2020). In our study, nine patients had AIS, and of these, only one patient, P7 who was a 31-year-old female, was diagnosed with gonadal malignancy. While evaluating the clinical details of the other patients, malignancy was not seen, which could be associated with younger age at the time of presentation; these patients need to be followed up to determine their status of malignancy.

All the patients with SRD5A2 gene mutations (Patients P10–P14) had a male karyotype, ambiguous genitalia, and variable degrees of under-virilization ranging from pseudo-vaginal hypospadias to clitoromegaly and variable positions of either one or both the gonads. Müllerian structures were found to be absent in all the patients. Although the ratios of testosterone to DHT were not available for the exact clinical diagnosis of these cases due to insufficient data, NGS was found to be important in confirming the diagnosis of steroid 5 alpha-reductase deficiency.

As far as the SF1 gene is concerned, all three cases showed a 46,XY karyotype; two of these cases, P15 and P17, had ambiguous genitalia, and one case, P16, had female external genitalia. The secondary sexual characters were typically underdeveloped with variable Tanner staging of development. The uterus was found to be absent in all three cases and either or both the ovaries were found to be small or absent in Patients P15 and P16. In Patient P17, a gonadectomy revealed the gonads to be testicular tissue. All the three patients with SF1 gene mutations had elevated levels of FSH and were categorized as hypergonadotropic hypogonadism.

A clinical diagnosis was difficult in Patients P18, P19, and P20, as all of them had a male karyotype, ambiguous genitalia, and absent uterus and ovaries with elevated levels of FSH. These patients were identified with mutations in the DHX37, DMRT2, and HSD17B3 genes, respectively, providing them with an exact molecular diagnosis, which will help in genetic counseling (Berra et al., 2011). Gonadal tumors occur in up to 25% of women with a Y chromosome; unlike cases of complete AIS, their gonads do not secrete hormones and should therefore be removed at the time of diagnosis. Although the management of patients with DSD is difficult due to the complex molecular mechanisms, the evaluation of tumor risk can be aided by the advances in genotyping for Y-chromosomal material, which is not evident in traditional karyotyping. Thus, the benefits of reaching a specific diagnosis are not only important for understanding the tumor risk but for identifying associated features, inheritance, chances of other family members being affected, and an understanding of the natural history of a condition. In our cohort of 25 patients, we could not find any mutation in the DSD-related genes in five patients (20%, 5/25) due to complex molecular mechanisms in 46,XY DSD; this suggests new DSD genes that are yet to be discovered in these disorders. Even though there are already known genes in DSD patients, it is difficult to establish the genotype–phenotype co-relation (Wang et al., 2017). Furthermore, the pathogenicity of novel variants was determined by multiple variant effect prediction tools and protein modeling and docking studies; however, in vitro functional studies are important to ascertain the pathogenicity of these variants, which can be considered as the limitation of this study.

Conclusion

The present study concludes that adolescent females with DSD and presenting with PA and absent secondary sexual characters should be investigated for chromosomal abnormalities along with routine hormonal and ultrasound investigations. Cytogenetically confirmed 46,XY adolescent females should be investigated for genetic mutations at a molecular level to arrive at a specific diagnosis, as the clinical signs and symptoms are found to be variable in these cases. Targeted NGS using a gene panel related to sex development will help in identifying an exact molecular diagnosis in these patients and should be preferred over routine Sanger sequencing. Our study highlights the importance of molecular diagnosis and genetic counseling of DSD patients, which will ultimately help in appropriate gender assignment of these cases.

Supplementary data

Supplementary data are available at Molecular Human Reproduction online.

Data availability

All data generated or analyzed during this study are included in this article and its supplementary information files. The novel variants identified have been reported to ClinVar in a two-spreadsheet submission (Accession no: SCV001960987–SCV001961002 and SCV002097635–SCV002097638).

Acknowledgments

We are grateful to ICMR-NIIH for providing financial support. We thank all the patients for participating in our study and the technical staff of Department of Cytogenetics, ICMR-NIIH.

Authors’ roles

V.K., S.K.C., S.D., and B.R.V. contributed to the study conception and design, carried out the experiments, and analyzed the data. V.K. performed the clinical assessment, diagnosis of patients, and experimental work. S.K.C. performed the bioinformatic analysis and in silico validation. V.K., S.D., and J.G. collected and analyzed the data and carried out the validation work. V.K. and S.K.C. prepared the article draft. B.R.V. edited and finalized the article. All the authors read and approved the article before submission.

Funding

This study was funded by the Indian Council of Medical Research through an Intramural grant provided to ICMR-NIIH (Intramural/ICMR-NIIH/2018-19).

Conflict of interest

The authors declare that there are no conflicts of interest in connection with this article.

References