A Boolean network model of human gonadal sex determination
© Ríos et al. 2015
Received: 21 August 2015
Accepted: 30 October 2015
Published: 16 November 2015
Gonadal sex determination (GSD) in humans is a complex biological process that takes place in early stages of embryonic development when the bipotential gonadal primordium (BGP) differentiates towards testes or ovaries. This decision is directed by one of two distinct pathways embedded in a GSD network activated in a population of coelomic epithelial cells, the Sertoli progenitor cells (SPC) and the granulosa progenitor cells (GPC). In males, the pathway is activated when the Sex-Determining Region Y (SRY) gene starts to be expressed, whereas in females the WNT4/ β-catenin pathway promotes the differentiation of the GPCs towards ovaries. The interactions and dynamics of the elements that constitute the GSD network are poorly understood, thus our group is interested in inferring the general architecture of this network as well as modeling the dynamic behavior of a set of genes associated to this process under wild-type and mutant conditions.
We reconstructed the regulatory network of GSD with a set of genes directly associated with the process of differentiation from SPC and GPC towards Sertoli and granulosa cells, respectively. These genes are experimentally well-characterized and the effects of their deficiency have been clinically reported. We modeled this GSD network as a synchronous Boolean network model (BNM) and characterized its attractors under wild-type and mutant conditions.
Three attractors with a clear biological meaning were found; one of them corresponding to the currently known gene expression pattern of Sertoli cells, the second correlating to the granulosa cells and, the third resembling a disgenetic gonad.
The BNM of GSD that we present summarizes the experimental data on the pathways for Sertoli and granulosa establishment and sheds light on the overall behavior of a population of cells that differentiate within the developing gonad. With this model we propose a set of regulatory interactions needed to activate either the SRY or the WNT4/ β-catenin pathway as well as their downstream targets, which are critical for further sex differentiation. In addition, we observed a pattern of altered regulatory interactions and their dynamics that lead to some disorders of sex development (DSD).
KeywordsSex determination Gonadal sex determination Boolean model Gene regulatory network
Sex development is a complex biological process that occurs during the embryonic and fetal stages of an individual. For a better understanding sex development is divided into three consecutive steps: 1) chromosomal sex determination (CSD); 2) gonadal sex determination (GSD); and 3) phenotypic sex differentiation (PSD). CSD is established at conception when the complement of sex chromosomes, XX or XY, is received. GSD, which is the process that we analyze in this study, refers to the set of genes and their regulatory interactions that trigger the development toward testes or ovaries, underlined by a gene regulatory network [1–4]. Finally, PSD involves the development of the female and male internal and external genitalia in response to the hormones secreted by the ovaries and testes. Both male and female PSD occur in two temporal phases, the first occurs within the fetus after GSD and the second occurs during puberty [5, 6].
GSD occurs within a heterogeneously composed structure called bipotential gonadal primordium (BGP). This structure, located on the ventromedial surface of the mesonephros [5–7], is critical for sex development since it can differentiate either as testes or ovaries . The BGP originates the actual gonad that is composed by a) the germinal cells (GCs), b) the steroidogenic somatic cells, such as the theca cells in ovary and the Leydig cells in testis that produce stradiol and testosterone, respectively; and c) the support somatic cells, including granulosa cells in ovary and Sertoli cells in testis.
Sertoli and granulosa cells originate from a common population of coelomic epithelial cells corresponding to the Sertoli progenitor cells (SPC) or granulosa progenitor cells (GPC) that migrate towards the BGP [9, 10]. In males, the SPCs start to differentiate toward Sertoli cells after 44 days of development (Carnegie-Stage 18). The mechanism involves activation of the expression of the Sex-determining Region Y gene (SRY) that codifies the SRY transcription factor [9, 11]. SRY associates with other transcription factors (i.e., CBX2, SF1) to regulate expression of the SOX9 gene that positively regulates the expression of genes associated to Sertoli cells (i.e., SOX9, FGF9, PGD2, DHH, AMH) . In females, where SRY is absent, GSD initiates after 49 days of development (Carnegie-Stage 20). In this case, the GPCs of coelomic origin differentiate towards granulosa cells by the action of a distinct gene regulatory pathway. Most likely, an increased amount of the transcription factor β-catenin up-regulates a set of downstream genes associated to granulosa, such as FOXL2 and RSPO1 [2, 12, 13]. Thus, the mechanism underlying GSD involves a common population of undifferentiated cells with the potential to diverge into two cell fates. The male pathway leads towards Sertoli cell fate determination and differentiation, whereas the female pathway leads to granulosa cell fate determination and differentiation.
Once differentiated, the Sertoli cells act as organizing centers, enclosing GCs to form testicular cords and secreting factors such as DHH and PDGF, which are essential for development of the fetal population of Leydig cells . Granulosa cells are the female equivalent of the Sertoli cells, as they enclose GCs and secrete factors necessary for oocyte growth and maturation. The regulatory network controlling GSD and differentiation toward Sertoli or granulosa cell consists, in a broad sense, of multiple target genes, different types of RNAs, transcription factors, nuclear receptors and signaling molecules. These elements are present in undifferentiated cells and interact in a concerted way either activating or repressing target genes at the time of GSD to balance the fate toward Sertoli or granulosa cells [15–17].
The total number of genes implicated in the regulatory network of GSD of humans and mammals remains elusive, as well as their complete regulatory interactions and their effects on the process of Sertoli or granulosa cells differentiation. However, it is well known that mutations in their components underlay the so-called disorders of sex development (DSD), a series of genetic conditions characterized by anomalies in gonads as well as in internal and external genitalia. The incidence of DSDs, as estimated by the The Lawson Wilkins Pediatric Endocrine Society (LWPES) and the European Society for Pediatric Endocrinology (ESPE), is 1 in 4,500 births  and can be attributed to mutations in various genes of the GSD network. For example, mutations in CBX2, GATA4 and WT1 genes result in a wide range of phenotypic alterations characterized by ambiguous or female external genitalia with the presence or absence of Mullerian structures in 46,XY DSDs patients . In contrast, 46,XX DSDs cause masculinization of the female fetus (normal males with no ovarian tissue) . In other cases 46,XX DSD patients have a female phenotype but fail to develop ovaries, presenting instead a “streak gonad”? (streaks of connective fibrous tissue) .
Boolean network models (BNM) are formal tools for analyzing the structure and dynamic behavior of genetic regulatory networks. BNM are best suited for describing poorly-characterized systems with no or few kinetic details, such as the GSD network. These models represent molecular entities (genes, transcription factors and RNAs) as nodes interacting among them within a network. Each node can have only two qualitative states: 0 (OFF) and 1 (ON). The OFF state is equivalent to a below-threshold concentration or activity, which is insufficient to initiate the intended process or regulation, while the ON state is equivalent to an above-threshold concentration or activity . The ON/OFF state of a node within the network is determined by a Boolean function that encompasses the known regulatory elements of the target node (transcription factors, nuclear receptors, signaling molecules). The state of these regulatory elements is updated over consecutive time steps of a simulation until the system converges to either a steady state or a cycle. BNMs describe the dynamic state of the nodes in a network by updating the state of the nodes according with the set of regulatory functions . BNMs have been implemented for the analysis of developmental programs such as flower morphogenesis in A. thaliana , early cardiac development in mice , and expression pattern of the segment polarity genes in Drosophila  to name a few.
Despite the relatively high incidence of DSDs, their molecular basis at the level of the regulatory network remain poorly understood. Thus, we are interested in constructing a BNM of the process of gonadal sex determination with an emphasis on the regulatory elements that are present at early stages of development and control the differentiation of SCP and GCP towards Sertoli and granulosa cells, respectively, allowing us to analyze the origin of some DSDs. For in-depth reviews about the genes involved in GSD and DSDs see: [18, 19, 27, 28], as well as the list of genes and interactions in the Additional file 1 of the supplementary information.
In this study we present a BNM that describes the dynamics of the GSD regulatory network starting from the UGR until Sertoli/granulosa cells differentiation. The proposed regulatory network incorporates a large amount of published information related to functional interactions among the genes involved in this process, while the BNM of GSD describes the dynamics of the elements contained within the network under wild-type and mutant conditions. With the current model we explore a formal description of the functional relationships among the genes and gene products associated to GSD, and generate some predictions about the expected regulatory behavior under wild-type or altered conditions within elements of the UGR and elements of the bipotential gonadal primordium such as CBX2, GATA4, and WT1. Additional predictions are indicated in the female pathway where the transcription factor β-catenin seems to play an important role in the activation of female-specific genes (for example, WNT4, RSPO1, and FOXL2).
The network of gonadal sex determination
The regulatory interactions among nodes were inferred, with emphasis in humans, from: (1) clinical studies of patients with DSDs, which carried mutations on sex determining genes; (2) genetic expression patterns associated to GSD (between 5th and 8th weeks of embryonic development); and (3) molecular evidence of interactions at the level of transcriptional regulation of target genes (i.e., up or down regulation of a target gene by means of protein-DNA interactions under wild type, mutated or transgenic constructs). Experimental evidence on mice was integrated into the network when necessary, especially in the female pathway where human information is lacking. References to clinical and experimental data can be found in the Additional file 1 of the Supplementary Information.
The network includes the special UGR node, representing the urogenital ridge, an embryonic structure precursor of the nephrogenic cord and gonads. The UGR node encompasses the following genes: LHX1, LHX9, EMX2, PAX2 and, PAX8. Although expression of these genes is essential for growth and maintenance of the UGR, little evidence was found in human, as well as in mouse, about their specific regulatory interactions, thus these genes were grouped within the UGR node, since mutations in any of these genes impair subsequent gonadal development [29–31].
The pathway towards pre-Sertoli or pre-granulosa cells is shown in Fig. 1. The blue nodes correspond to Sertoli cell fate determination pathway and include: SRY, SOX9, FGF9, PGD2, DHH, AMH, DKK1 and DMRT1 nodes, whereas the pink nodes correspond to granulosa cell fate determination pathway including: CTNNB1, WNT4, FOXL2 and RSPO1. Notice that the granulosa cell fate determination was complemented with mice information. For example, we considered the canonical Wnt4/ β-catenin pathway as a key regulatory element of female nodes within the network since relative expression of Fst, Gng13, Foxl2, Irx3 and, Sp5 has been shown to be down-regulated when β-catenin is lost in female mice in early stages of ovarian development .
The network as a Boolean model
The process of GSD is poorly characterized at the quantitative level, i.e., kinetic information regarding the interactions of the elements of this regulatory network is still lacking, therefore the implementation of the GSD network as a continuous model is, at this moment, out of reach. Given this, we decided to model the network as a discrete dynamical system so as to describe the qualitative observations that are experimentally reported. Specifically, we used a Boolean approach where every node might have one of two possible states; 1 (ON) or 0 (OFF), indicating that a given node within the network model is active or inactive, respectively.
Set of functions for the Boolean model of gonadal sex determination
UGR, UGR & ! (NR5A1 ∣ WNT4)
CBX2, UGR & ! (NR0B1 & WNT4 & CTNNB1)
GATA4, (UGR ∣ WNT4 ∣ NR5A1 ∣ SRY)
WT1mKTS, (UGR ∣ GATA4)
WT1pKTS, (UGR ∣ GATA4) & ! (WNT4 & CTNNB1)
NR5A1, (UGR ∣ CBX2 ∣ WT1mKTS ∣ GATA4) & ! (NR0B1 & WNT4)
NR0B1, (WT1mKTS ∣ (WNT4 & CTNNB1)) & ! (NR5A1 & SOX9)
SRY, ((NR5A1 & WT1mKTS & CBX2) ∣ (GATA4 & WT1pKTS & CBX2 & NR5A1) ∣ (SOX9 ∣ SRY)) & ! (CTNNB1)
SOX9, ((SOX9 & FGF9) ∣ (SRY ∣ PGD2) ∣ (SRY & CBX2) ∣ (GATA4 & NR5A1 & SRY)) & ! (WNT4 ∣ CTNNB1 ∣ FOXL2)
FGF9, SOX9 & ! WNT4
DMRT1, (SRY ∣ SOX9) & ! (FOXL2)
DKK1, (SRY ∣ SOX9)
AMH, ((SOX9 & GATA4 & NR5A1) ∣ (SOX9 & NR5A1 & GATA4 & WT1mKTS)) & ! (NR0B1 & CTNNB1)
WNT4, (GATA4 ∣ (CTNNB1 ∣ RSPO1 ∣ NR0B1)) & ! (FGF9 ∣ DKK1)
RSPO1, (WNT4 ∣ CTNNB1) & ! (DKK1)
FOXL2, (WNT4 & CTNNB1) & ! (DMRT1 ∣ SOX9)
CTNNB1, (WNT4 ∣ RSPO1) & ! (SRY ∣ (SOX9 & AMH))
We performed an initial exhaustive evaluation of the dynamic behavior of the wild type model, simulating all possible initial activation states. Three fixed-point attractors were obtained, and we performed a search focused in finding the state transitions corresponding to both male and female pathways. To recover the wild type “male pathway”, we initiated the simulations with the UGR node in ON. In contrast, to created a wild type “female pathway”, without the SRY node, we set the UGR and WNT4 nodes as active at the beginning of simulations. Besides the wild type model, we simulated all possible loss and gain of function of single mutants, so as to describe alterations in activation states that might be interpreted as alterations in gene expression. Loss and gain of function single mutants were simulated by fixing the relevant node to 0 or 1, respectively. All simulations were carried out under the synchronous updating scheme with the use of BoolNet .
Testing properties of the Boolean model: random networks and robustness of attractors
The network of gonadal sex determination
The network was constructed with 19 nodes and 78 regulatory interactions: 42 of these interactions have been reported in humans; 12 in mice and 24 were predicted from analysis of transition states of the simulated Boolean model. The network is directed towards the male pathway if the SRY node is active. SRY leads to activation of SOX9, which in turn activates FGF9, PGD2, DMRT1, DHH, DKK1, and AMH nodes. At the same time, the female pathway is repressed by inactivating CTNNB1 and FOXL2 nodes (Fig. 1). On the contrary, the network is directed towards the female pathway in absence of SRY and when the WNT4, CTNNB1, RSPO1 and, FOXL2 nodes are active. In this case, the male pathway is repressed by CTNNB1 and FOXL2 mediated inactivation of SOX9, DMRT1 and AMH (Fig. 1).
Predictions of the Boolean model
The current model contains 24 interactions inferred from dynamic modeling, these are predominantly related to the UGR node and the genes expressed in the bipotential gonad. Model predictions were drawn as orange dashed lines within the following nodes: UGR, CBX2, GATA, WT1mKTS, WT1pKTS, NR0B1, SRY, DMRT1, DKK1, WNT4 and CTNNB1 (Fig. 1). The model predicts that the activity of UGR depends of an activation self-loop and functions as an input to activate CBX2, GATA4, Wt1mKTS, WT1pKTS, and NR5A1 (Table 2). The scarcity of information regarding UGR function and maintenance clearly indicates that more experimental studies are necessary to understand the mechanisms of gene expression control in the BGP especially for CBX2, GATA4 and WT1.
Dynamic behavior of the gonadal sex determination Boolean network model
The dynamic behavior of the GSD BNM was exhaustively analyzed by starting the dynamical simulations of the system from all possible 219=524288 initial states. After simulations, three fixed-point attractors where obtained (Fig. 2). The first of these attractors can be interpreted as the gene expression profile observed in Sertoli cells, the second can be interpreted as the gene expression profile observed in granulosa cells, and the third attractor, with a very small basin of attraction, might represent a disgenetic gonad without Sertoli or granulosa activity.
If the simulation transited toward the Sertoli attractor, then the NR0B1, WNT4 and RSPO1 nodes were inactivated by NR5A1, SRY, SOX9 and DKK1 nodes. In male the expression of NROB1 is dosage sensitive, since duplication of this gene in 46,XY patients produces a male-to-female sex reversal with streak gonads . It is important to notice that NR0B1 might play an important role in male after the time of GSD because NR0B1 knockout mice showed disorganized Sertoli, Leydig and germ cells due to defects in testis cord formation . Thus, it has been suggested that NR0B1 has a time frame of expression  with reduced levels of the DAX1 protein during GSD. Since the BNM considers active or inactive states, the NR0B1 node was inactive at the sixth time step, which corresponds to the Sertoli attractor (Fig. 3 a)
When the UGR + WNT4 nodes were set to ON, two fixed point attractors were obtained: (1) the granulosa attractor and the (2) dysgenetic gonad attractor. Thus the transition towards the granulosa attractor is characterized by the initial activation of the UGR + WNT4 and BPG nodes, followed by activation of the NR0B1, RSPO1, FOXL2 and CTNNB1. As we previously stated, β-catenin, plays a key role in up-regulation of pre-granulosa genes in female mouse , this factor actively antagonizes SOX9 and AMH expression, inactivating the pathway toward Sertoli cells (Fig. 3 b) [38, 39]. The attractor with no activity reflects the importance of the UGR node within the network model given that loss of function mutants of UGR components have an impaired subsequent gonadal development, as observed in mouse. Therefore, the dysgenetic gonad attractor (i.e., streaks of fibrous tissue instead of a gonad) might be interpreted as a condition expected in some individuals when gonadal development fails, especially in the case of LHX1, LHX9, EMX2, PAX2 and PAX8 mutants.
In summary, state transitions in Fig. 3 a and b qualitatively coincide with gene expression patterns observed during GSD [7, 34]. However, notice that the state transitions and steady state attractors must be considered as snapshots of the gene expression pattern between 41–52 days of development and do not represent the complete process of gonadal development.
Modeling disorders of sex development
46,XX sex reversal
46,XY (SRY-) sex reversal
SRY is considered the trigger of the male pathway and testis development . This transcription factor possesses a highly conserved domain called HMG box that binds to the GAACAAAG DNA motif and bends the DNA molecule about 80 degrees. The loss of its chromatin-remodeling activity  is considered to impair the three dimensional architecture of chromatin and compromises the proper interaction of SRY with its target genes. Mutation of the DNA binding region of SRY in 46,XY subjects has been associated with female external genitalia, normal Mullerian ducts and streak gonads . When we simulated the SRY loss of function, the female pathway was activated by the CTNNB1 node and the male pathway blocked through a CTNNB1 and FOXL2- mediated SOX9 inhibition. SRY is considered the trigger of testis development by expressing in the somatic pre-Sertoli cells . The mechanism suggested for normal function of this transcription factor is a highly conserved domain within the protein, called High Mobility Group (HMG box).
Modeling other DSDs
We analyzed in addition the effect of loss-of-function mutations in the NR5A1 node. Since the dynamics of the simulations of WT1 and NR5A1 were identical we show only the attractor in Fig. 7.
Testing properties of the Boolean model: random networks and robustness of the attractors
On average, the three independent tests of 1000 random networks generated 103,000 attractors each one. We found that none of the random networks recovered the set of three point attractors of the wild type model shown in Fig. 2. This results indicate that the attractors of our BNM could not be expected in a network made with random interactions. Thus, the attractors in Fig. 2 can be considered as biologically meaningful and not a statistical artifact.
Concerning the test of attractor robustness we observed that 95 % of the perturbed networks recovered less than 30 % of the original attractors of our model. This means that the model is relatively sensitive to perturbations in the functions of the model. The reduced robustness can be due to the scarce redundancy in the model, given that we opted by including a small number of regulatory molecules, so as to be close to a minimal model. We expect that the introduction of more nodes and regulatory feedback circuits would result in an increased robustness.
The classical observations of Alfred Jost (1947) that early castration in utero of rabbit fetuses resulted in female internal and external genitalia (independent of their chromosomal sex complement) lead to the hypothesis of a testis-determining factor (TDF). According to this hypothesis the ovary was considered the default developmental state, while testis represented an induced and active state that repressed female development. Experimental evidence accumulated in the last 20 years have enriched our view of sex determination where developmental programs toward testis or ovaries represent two independent antagonistic regulatory pathways of high complexity intertwined in a regulatory network: the GSD network.
The BNM used in this study, although discrete in its approach, can be considered as a simplified version of a very dynamic and complex biological network that incorporates the major regulatory elements of the GSD network. The GSD BNM summarizes in a formal language the set of experimentally-confirmed interactions associated with the process of GSD. Although the attractors obtained in our model cannot be interpreted as as anatomical structures of high developmental complexity, they can be reliably seen as the gene expression profiles expected during the process of determination and differentiation of SPC and GPC towards Sertoli and granulosa cells respectively.
The network of gonadal sex determination
We inferred the regulatory network of human GSD and modeled it as a BNM. With such a model we were able to describe the molecular dynamics of the first stage in gonadal morphogenesis, which is the cell fate determination and further differentiation of Sertoli and granulosa cells. The network contains 19 nodes as well as their regulatory interactions, as evidenced by published experimental and clinical data.
Recent studies regarding early gonadal differentiation suggest a highly complex biological process regulated by many, probably hundreds, of genes. However, our BNM shows that only a handful of them are sufficient to activate the male or female pathway, allowing us to propose that the gonadal fate commitment and differentiation is a direct consequence of activation and repression of a transcriptional program encoded as a regulatory network. Although part of the information used to infer this regulatory network was taken from experiments in mice, we consider that the model might be considered as a good approximation to the corresponding regulatory network in humans.
The Boolean model of gonadal sex determination
The BNM describes the dynamics of 19 nodes associated with GSD between 41–52 days of embryonic development. The results cannot be considered as final activation states of the biological process, instead they should be considered as a snapshot in the process of Sertoli and granulosa cell differentiation. At the quantitative level, GSD is poorly understood given the lack of information about kinetic details of each regulatory element, therefore it is difficult to establish a continuous model with differential equations with the current available data. Despite the large amount of gene expression data, little is known about the regulatory mechanisms leading to GSD under normal conditions, as well as their downstream effects under mutant conditions. To shed light about the regulatory mechanisms we used a discrete modeling approach because most of the information relies on qualitative descriptions.
Model predictions and the role of the UGR and bipotential gonad genes in early gonadal development
The UGR is a very important structure within the developing embryo since it is the common structure that leads to testis and ovaries. Notice that just a few regulatory elements of the UGR have been studied. For example, Lhx1 expression has been reported in mice and has a key role in the development of kidney, female reproductive tract and anterior head . Conditional knockout mice lack uterus, cervix and upper vagina . Lhx9 has a role in the activation of Nr5a1 gene, in synergy with Wt1 in mice. Other example is given by PAX2 and PAX8 genes. In humans, these genes have a role in the activation of the WT1 gene [29, 50] thus, we underline the need for additional studies regarding the regulatory interactions that led to the establishment of the UGR and BPG primordium.
As we mentioned previously, the human female pathway is less characterized, and its current cumulative knowledge is mainly based on mice findings. In this case, the GATA4-FOG2 complex has an important function by activating the Fst, Wnt4, Sprr2d, Foxl2, Gng13 genes . It has been observed in female mice that loss of function of Gata4 impairs the expression of these genes and leads to the development of a male-specific coelomic vessel ; therefore GATA4 can be considered as a key regulatory element in the early stages of gonadal development toward ovaries.
The role of WNT4/ β-catenin in the female pathway
Concerning the female pathway most of the inferred interactions derive from mice. To this respect we notice the role of the WNT4 and βcatenin nodes in the regulation of WNT4, RSPO1, and FOXL2. In the biological process we predict a key role of β-catenin regulating the female pathway. The general mechanism of the canonical Wnt4/ β-catenin signaling pathway can be explained as follows: the pathway is initiated by expression of the wingless-type MMTV integration site family, member 4 (WNT4). The product of this gene is a ligand that binds to Frizzled (Fz) and LRP5/6 co-receptors at the plasma membrane, disengaging β-catenin from the proteins of the “destruction complex”? (Axin and APC). Then β-catenin translocate into the nucleus where associates with TCF7/LEF, this protein contains an HMG box with capacity to recognize specific DNA sequences. The β-catenin-TCF7/LEF complex activates target genes, whereas in absence of β-catenin TCF7/LEF alone represses gene transcription [51, 52]. From the set of interactions inferred from mice, that fitted perfectly in our model, we would expect that β-catenin might have an important role regulating the expression of the FST, FOXL2 and IRX3 genes in humans.
The BNM of GSD summarizes in a formal language the set of experimentally-confirmed interactions associated with the process of GSD. The attractors of our model can be interpreted as the gene expression profiles expected during the process of GSD and differentiation of Sertoli or granulosa cell lineages. According to our simulations, the loss of function of GATA4 results in inactivation of the AMH node in the attractor. This result is particularly interesting given the existence of the persistent Müllerian duct syndrome (PMDS), a relatively rare inherited defect in the sexual differentiation, characterized by failure in the regression of the Müllerian ducts in males. Affected individuals present persistent uterus and tubes due to a defect in the synthesis of the AMH hormone, which is normally produced by the Sertoli cells. Mutations in the AMH gene have been reported in these patients [20, 53], however the majority of PMDS remain without molecular diagnosis, therefore GATA4 mutations emerge, according to our model predictions, as a potential PMSD causing gene.
We propose the present model as a starting point for future mathematical modeling and integration of experimental research regarding sex development. The model can be upgraded in several aspects for example, incorporating additional nodes and interactions, as well as modeling more cell lineages of the gonad such as the Leydig or theca cells. Finally the current BNM describes the dynamics of the GSD network under perturbations. Importantly the analysis of these states can have potential implications in the study of DSDs.
OR thanks the support from the program Doctorado en Ciencias Biológicas, UNAM, and the scholarship 173000 from Consejo Nacional de Ciencia y Tecnología (CONACyT). This work was supported by grants from CONACYT 166012 to HM, Universidad Nacional Autónoma de México PAPIIT IN200514 to LM, and Recursos Fiscales para la Investigación del Instituto Nacional de Pediatría, proyecto 057/2014 to LT. Adhemar Liquitaya from CompBioLab IIB, UNAM, aid with the test of random networks.
Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http://creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated.
- De Santa Barbara P, Moniot B, Poulat F, Berta P. Expression and subcellular localization of SF-1, SOX9, WT1, and AMH proteins during early human testicular development. Dev Dynam. 2000; 217(3):293–8.View ArticleGoogle Scholar
- Eggers S, Sinclair A. Mammalian sex determination insights from humans and mice. Chromosome Res. 2012; 20(1):215–38.PubMed CentralView ArticlePubMedGoogle Scholar
- White S, Sinclair A. The molecular basis of gonadal development and disorders of sex development In: Hutson JM, Warne GL, Grover SR, editors. Disorders of sex development: An Integrated Approach to Management. Heidelberg: Springer Berlin: 2012. p. 1–9.Google Scholar
- She ZY, Yang WX. Molecular mechanisms involved in mammalian primary sex determination. J Mol Endocrinol. 2014; 53(1):R21–37.View ArticlePubMedGoogle Scholar
- Sajjad Y. Development of the genital ducts and external genitalia in the early human embryo. J Obstet Gynaecol Re. 2010; 36(5):929–37.View ArticleGoogle Scholar
- Hutson JM. Embryology of the human genital tract In: Hutson JM, Warne GL, Grover SR, editors. Disorders of sex development: An Integrated Approach to Management. Heidelberg: Springer Berlin: 2012. p. 11–21.View ArticleGoogle Scholar
- Hanley NA, Hagan DM, Clement-Jones M, Ball SG, Strachan T, Salas-Cortes L, et al.SRY, SOX9, DAX1 expression patterns during human sex determination and gonadal development. Mech Develop. 2000; 91(1):403–7.View ArticleGoogle Scholar
- Kim Y, Capel B. Balancing the bipotential gonad between alternative organ fates: a new perspective on an old problem. Dev Dynam. 2006; 235(9):2292–300.View ArticleGoogle Scholar
- Albrecht KH, Eicher EM. Evidence that SRY is expressed in pre-Sertoli cells and Sertoli and granulosa cells have a common precursor. Dev Biol. 2001; 240(1):92–107.View ArticlePubMedGoogle Scholar
- Ungewitter EK, Yao HC. How to make a gonad: cellular mechanisms governing formation of the testes and ovaries. Sex Dev. 20013; 7(1-3):7–20.View ArticlePubMedGoogle Scholar
- Kashimada K, Koopman P. SRY: the master switch in mammalian sex determination. Development. 2010; 137(23):3921–30.View ArticlePubMedGoogle Scholar
- Harikae K, Miura K, Kanai Y. Early gonadogenesis in mammals: significance of long and narrow gonadal structure. Dev Dynam. 2013; 242(4):330–8.View ArticleGoogle Scholar
- Sarraj MA, Drummond AE. Mammalian foetal ovarian development: consequences for health and disease. Reproduction. 2012; 143(2):151–63.View ArticlePubMedGoogle Scholar
- Svechnikov K, Landreh L, Weisser J, Izzo G, Colon E, Svechnikova I, et al.Origin development and regulation of human Leydig cells. Horm Res Paediatr. 2010; 73(2):93–101.View ArticlePubMedGoogle Scholar
- Beverdam A, Koopman P. Expression profiling of purified mouse gonadal somatic cells during the critical time window of sex determination reveals novel candidate genes for human sexual dysgenesis syndromes. Hum Mol Genet. 2006; 15(3):417–31.View ArticlePubMedGoogle Scholar
- Jameson SA, Natarajan A, Cool J, DeFalco T, Maatouk DM, Mork L, et al.Temporal transcriptional profiling of somatic and germ cells reveals biased lineage priming of sexual fate in the fetal mouse gonad. PLoS Genet. 2012; 8(3):e1002575.PubMed CentralView ArticlePubMedGoogle Scholar
- Munger SC, Natarajan A, Looger LL, Ohler U, Capel B. Fine time course expression analysis identifies cascades of activation and repression and maps a putative regulator of mammalian sex determination. PLoS Genet. 2013; 9(7):e1003630.PubMed CentralView ArticlePubMedGoogle Scholar
- Hughes IA, Houk Ch, Ahmed SF, Lee PA. Consensus statement on management of intersex disorders. J Pediatr Urol. 2006; 2(3):148–62.View ArticlePubMedGoogle Scholar
- Mendonca BB, Domenice S, Arnhold IJ, Costa EM. 46,XY disorders of sex development (DSD). Clin Endocrinol. 2009; 70(2):173–87.Google Scholar
- López M, Torres L, Méndez JP, Cervantes A, Perez-Palacios G, Erickson RP. Clinical traits and molecular finding in 46,XX males. Clin Genet. 1995; 48(1):29–34.View ArticlePubMedGoogle Scholar
- Kousta E, Papathanasiou A, Skordis N. Sex determination and disorders of sex development according to the revised nomenclature and classification in 46, XX individuals. Hormones (Athens). 2010; 9(3):218–31.View ArticleGoogle Scholar
- Saadatpour A, Albert R. Boolean modeling of biological regulatory networks: a methodology tutorial. Methods. 2013; 62(1):3–12.View ArticlePubMedGoogle Scholar
- Karlebach G, Shamir R. Modelling, analysis of gene regulatory networks. Nat Rev Mol Cell Bio. 2008; 9(10):770–80.View ArticleGoogle Scholar
- Sanchez-Corrales YE, Alvarez-Buylla ER, Mendoza L. The Arabidopsis thaliana flower organ specification gene regulatory network determines a robust differentiation process. J Theor Biol. 2010; 264(3):971–83.View ArticlePubMedGoogle Scholar
- Herrmann F, Groß A, Zhou D, Kestler HA, Kühl M. A Boolean model of the cardiac gene regulatory network determining first and second heart field identity. PLoS ONE. 2012; 7(10):e46798.PubMed CentralView ArticlePubMedGoogle Scholar
- Albert R, Othmer HG. The topology of the regulatory interactions predicts the expression pattern of the segment polarity genes in Drosophila melanogaster. J Theor Biol. 2003; 223(1):1–18.View ArticlePubMedGoogle Scholar
- Fleming A, Vilain E. The endless quest for sex determination genes. Clin genet. 2005; 67(1):15–25.View ArticlePubMedGoogle Scholar
- Ono M, Harley VR. Disorders of sex development: new genes new concepts. Nat Rev Endocrinol. 2013; 9(2):79–91.View ArticlePubMedGoogle Scholar
- McConnell MJ, Cunliffe HE, Chua LJ, Ward TA, Eccles MR. Differential regulation of the human Wilms tumour suppressor gene (WT1) promoter by two isoforms of PAX2. Oncogene. 1997; 14(22):2689–700.View ArticlePubMedGoogle Scholar
- Kobayashi A, Shawlot W, Kania A, Behringer RR. Requirement of Lim1 for female reproductive tract development. Development. 2004; 131(3):539–49.View ArticlePubMedGoogle Scholar
- Wilhelm D, Englert Ch. The Wilms tumor suppressor WT1 regulates early gonad development by activation of SF1. Genes Dev. 2002; 16(14):1839–51.PubMed CentralView ArticlePubMedGoogle Scholar
- Manuylov NL, Smagulova FO, Leach L, Tevosian SG. Ovarian development in mice requires the GATA4-FOG2 transcription complex. Development. 2008; 135(22):3731–43.View ArticlePubMedGoogle Scholar
- Müssel Ch, Hopfensitz M, Kestler HA. BoolNet an R package for generation reconstruction and analysis of Boolean networks. Bioinformatics. 2010; 26(10):1378–80.View ArticlePubMedGoogle Scholar
- Tomaselli S, Megiorni F, Lin L, Mazzilli MC, Gerrelli D, Majore S, et al.Human RSPO1/R-spondin1 is expressed during early ovary development and augments β-catenin signaling. PLoS One. 2011; 6(1):e16366.PubMed CentralView ArticlePubMedGoogle Scholar
- Barbaro M, Oscarson M, Schoumans J, Staaf J, Ivarsson SA, Wedell A. Isolated 46, XY gonadal dysgenesis in two sisters caused by a Xp21. 2 interstitial duplication containing the DAX1 gene. J Clin Endocrinol Metab. 2007; 92(8):3305–13.View ArticlePubMedGoogle Scholar
- Meeks JJ, Crawford SE, Russell TA, Morohashi K, Weiss J, Jameson JL. Dax1 regulates testis cord organization during gonadal differentiation. Development. 2003; 130(5):1029–36.View ArticlePubMedGoogle Scholar
- Ludbrook LM, Harley VR. Sex determination: a window of DAX1 activity. Trends Endocrin Met. 2004; 15(3):116–21.View ArticleGoogle Scholar
- Hersmus R, Kalfa N, de Leeuw BHCGM, Stoop H, Oosterhuis JW, de Krijger R, et al.FOXL2 and SOX9 as parameters of female and male gonadal differentiation in patients with various forms of disorders of sex development (DSD). J Pathol. 2008; 215(1):31–8.View ArticlePubMedGoogle Scholar
- Maatouk DM, DiNapoli L, Alvers A, Parker KL, Taketo MM, Capel B. Stabilization of β-catenin in XY gonads causes male-to-female sex-reversal. Hum Mol Genet. 2008; 17(19):2949–55.PubMed CentralView ArticlePubMedGoogle Scholar
- Alves C, Braid Z, Coeli FB, Mello MP. 46 XX Male-testicular disorder of sexual differentiation (DSD): hormonal molecular and cytogenetic studies. Arq Bras Endocrinol. 2010; 54(8):685–9.View ArticleGoogle Scholar
- Harley VR, Clarkson MJ, Argentaro A. The molecular action and regulation of the testis-determining factors SRY (Sex-determining Region on the Y chromosome) and SOX9 [SRY-related high-mobility group (HMG) box 9]. Endocr rev. 2003; 24(4):466–87.View ArticlePubMedGoogle Scholar
- Hawkins JR, Taylor A, Goodfellow PN, Migeon CJ, Smith KD, Berkovitz GD. Evidence for increased prevalence of SRY mutations in XY females with complete rather than partial gonadal dysgenesis. Am J Hum Genet. 1992; 51(5):979.PubMed CentralPubMedGoogle Scholar
- Yueh-Chiang Hu, Okumura LM, Page DC. Gata4 is required for formation of the genital ridge in mice. PLoS Genet. 2013; 9(7):e1003629.View ArticleGoogle Scholar
- Miyamoto Y, Taniguchi H, Hamel F, Silversides DW, Viger RS. A GATA4/WT1 cooperation regulates transcription of genes required for mammalian sex determination and differentiation. BMC Mol Biol. 2008; 9(1):44.PubMed CentralView ArticlePubMedGoogle Scholar
- Tremblay JJ, Viger RS. A mutated form of Steroidogenic factor 1 (SF-1 G35E that causes sex reversal in humans fails to synergize with transcription factor GATA-4. J Biol Chem. 2003; 278(43):42637–42.View ArticlePubMedGoogle Scholar
- Lourenço D, Brauner R, Rybczyńska M, Nihoul-Fékété C, McElreavey K, Bashamboo A. Loss-of-function mutation in GATA4 causes anomalies of human testicular development. Proc Natl Acad Sci U S A. 2011; 108(4):1597–602.PubMed CentralView ArticlePubMedGoogle Scholar
- Köhler B, Biebermann H, Friedsam V, Gellermann J, Maier RF, Pohl M, et al.Analysis of the Wilms’ tumor suppressor gene (WT1) in patients 46, XY disorders of sex development. Clin Endocrinol Metab. 2011; 96(7):E1131–6.View ArticleGoogle Scholar
- Kim J, Prawitt D, Bardeesy N, Torban E, Vicaner C, Goodyer P, et al.The Wilms’ tumor suppressor gene WT1 product regulates DAX-1 gene expression during gonadal differentiation. Mol Cellular Biol. 1999; 19(3):2289–99.View ArticleGoogle Scholar
- Huang C-C, Orvis GD, Kwan KM, Behringer RR. Lhx1 is required in Müllerian duct epithelium for uterine development. Dev Biol. 2014; 389(2):124–36.PubMed CentralView ArticlePubMedGoogle Scholar
- Fraizer GC, Shimamura R, Zhang X, Saunders GF. PAX8 regulates human WT1 transcription through a novel DNA binding site. J Biol Chem. 1997; 272(49):30678–87.View ArticlePubMedGoogle Scholar
- Cadigan KM, Waterman ML. TCF/LEFs and Wnt signaling in the nucleus. Cold Spring Harbor Perspectives Biol. 2012; 4(11):a007906.View ArticleGoogle Scholar
- Tevosian SG, Manuylov NL. To β or not to β Canonical β-catenin signaling pathway and ovarian development. Dev Dynam. 2008; 237(12):3672–80.View ArticleGoogle Scholar
- Mazen I, Hamid A, El-Gammal M, Aref A, Amr K. AMH gene mutations in two Egyptian families with persistent Müllerian duct syndrome. Sex Dev. 2011; 5(6):277–80.View ArticlePubMedGoogle Scholar