Primary sex determination of placental mammals: a modelling study uncovers dynamical developmental constraints in the formation of Sertoli and granulosa cells

Background Primary sex determination in placental mammals is a very well studied developmental process. Here, we aim to investigate the currently established scenario and to assess its adequacy to fully recover the observed phenotypes, in the wild type and perturbed situations. Computational modelling allows clarifying network dynamics, elucidating crucial temporal constrains as well as interplay between core regulatory modules. Results Relying on a comprehensive revision of the literature, we define a logical model that integrates the current knowledge of the regulatory network controlling this developmental process. Our analysis indicates the necessity for some genes to operate at distinct functional thresholds and for specific developmental conditions to ensure the reproducibility of the sexual pathways followed by bi-potential gonads developing into either testes or ovaries. Our model thus allows studying the dynamics of wild type and mutant XX and XY gonads. Furthermore, the model analysis reveals that the gonad sexual fate results from the operation of two sub-networks associated respectively with an initiation and a maintenance phases. At the core of the process is the resolution of two connected feedback loops: the mutual inhibition of Sox9 and ß-catenin at the initiation phase, which in turn affects the mutual inhibition between Dmrt1 and Foxl2, at the maintenance phase. Three developmental signals related to the temporal activity of those sub-networks are required: a signal that determines Sry activation, marking the beginning of the initiation phase, and two further signals that define the transition from the initiation to the maintenance phases, by inhibiting the Wnt4 signalling pathway on the one hand, and by activating Foxl2 on the other hand. Conclusions Our model reproduces a wide range of experimental data reported for the development of wild type and mutant gonads. It also provides a formal support to crucial aspects of the gonad sexual development and predicts gonadal phenotypes for mutations not tested yet. Electronic supplementary material The online version of this article (doi:10.1186/s12918-016-0282-3) contains supplementary material, which is available to authorized users.


Basics
The model presented in this work was defined using the logical formalism (1,2). In this framework, the gene regulatory network is represented in terms of a directed graph. Each node (gene or regulatory product) is assigned a discrete variable with a maximal level, which defines the highest qualitative functional level taken by the regulatory node (this maximal level equals 1 in the simplest, Boolean case). Whenever distinct functional concentrations of a regulatory product need to be considered, multilevel variables are used. Each arc embodies a regulatory interaction and is assigned a threshold. This threshold defines the smallest functional level of the interaction source for which the interaction is operative. Logical parameters qualitatively describe the effects of each interaction or combination of interactions controlling the states of the network nodes.
Thus, each node is associated to as many parameters as the number of combinations of interactions acting on this node. Per default, these parameters are set to zero. Defining non zero parameters amounts to single out the combinations of interactions enabling the activation of the targeted node and to qualitatively recover all documented wild type and mutant gene expression patterns. Parameters for the 1-cell network of Fig. 1 are given in the Supplementary Table S2. A state of the network is represented by a vector, which encompasses the current (discrete) node levels. Given a state, one can determine which interactions are operative and the values of the logical parameters indicate the nodes called to change their levels (i.e. those nodes whose current levels differ from the values of the logical parameters). In general, for a given state, all possible elementary transitions (i.e. switching of a single node level, to a neighbouring integer level) are considered, thus leading to as many outgoing transitions as updating calls (asynchronous updating, see Section 1.3).
Depending on the structure of the regulatory graph and on the values of the logical parameters, a model leads to a large but finite number of dynamical pathways. These are in turn represented in the form of a State Transition Graph (STG), where nodes represent states and (directed) edges represent transitions between states.
The software tool GINsim (3) was used for the simulation and logical analysis of the gene regulatory network. The model files for the 1-cell and 2-cell networks are available as a Supplementary Files F1 and F2, in SBML format (4). They will be also provided in the model repository of the GINsim software web site (http://ginsim.org), where the tool is freely available.

Simulation of genetic perturbations
In the logical framework, simulation of genetic perturbations is straightforward. A lossof-function mutation of a given gene implies that this gene produces a non-functional product (or no product at all), which amounts to assign the value zero to the corresponding variable and parameters. In contrast, the ectopic expression of a gene implies that this gene is expressed in an unregulated manner beyond its normal spatiotemporal expression domain. This can be accomplished by forcing the corresponding variable to take higher values (for detail of this formal treatment of mutations, see (5).

Updating schemes and Hierarchical Transition Graphs
Model dynamics are defined by the specification of initial state(s) and of an updating scheme, which refers to how network component levels (i.e., model variables) are updated. In a synchronous update scheme, at each simulation step, all the components that are called to change their levels are simultaneously updated. This defines a deterministic behaviour in which each state has at most one successor. In contrast, the asynchronous update scheme is considered to be more realistic since it accounts for different activation and/or inactivation time of the model components. This update amounts to generate all possible trajectories from the specified initial conditions: e.g.
when two component levels are changed, these updates are done asynchronously, defining two transitions towards alternative successor states. The resulting STG may then include bifurcation states from which, depending on the choice of the successor state, the final stable state may differ; in other words, in the absence of any further constraint, two alternative phenotypes are reachable.
For our model of the gene network controlling primary sex determination, constraints on conflicting updates (transitions) needed to be considered to get the final state match the expected phenotype (testis or ovary). Such constraints are expressed in terms of priorities, i.e. temporal orders between component updates (6). To uncover these priorities, we relied on a compact representation of the dynamics, called Hierarchical transition graphs (HTG) (Supplementary Fig. S2).
HTG have been introduced as compact representations revealing the attractors and their basins of attraction, as well as transient oscillatory behaviours (7). To briefly describe the compaction leading to these graphs, we first recall that a Strongly Connected Component (SCC) is defined as a maximal strongly connected subgraph (i.e., a maximal subset of states in the STG, such that there is a path connecting each state to any other state). These SCCs are termed trivial when they encompass a single state, complex otherwise.
States of an STG are gathered into single nodes of the corresponding HTG as follows: • If they belong to a complex SCC; • If they belong to an attractor i.e. a trivial or a complex SCCs with no outgoing transitions; • If they define trivial SCCs from which the same complex and terminal SCCs are reachable.
Hence nodes of an HTG are defined as sets of states. The arcs between these nodes denote the existence of (at least) a transition between two states of the corresponding sets of states. They can be labelled by the associated updates of variables of the logical model (Supporting Fig. 2). Importantly, a path in a HTG towards any node defined by a complex or terminal SCC implies a path in the original STG (7).

Gene network controlling primary sex determination in placental mammals
In this section, the components and the interactions included in the gene network are discussed. The sexual cyto-differentiation genes, such as Amh and Fst were removed since they constitute the output, male versus female, of the regulatory network. Pgd2 was also eliminated because it is not essential for testis development (8). Gonad development into ovary requires the function of the canonical Wnt4 signalling pathway. This is triggered by both Wnt4 and Rspo1, which operate through the same effector, ß-catenin. For simplification, we selected the secreted ligand Wnt4 to represent Wnt4 signalling pathway.

Simplification of the regulatory network
The gene Sox9 plays a key, crucial role for testis determination, yet its function seems to be dispensable for maintaining adult testis identity (9). It appears that this function is exerted by its homolog Sox8, which is activated by Sox9 after its Sry-dependent upregulation (10). This is so because the elimination of Sox8 function in XY mice does not affect the development of the testis (11), except that they are unfertile, suggesting that Sox8 is required for the maintenance of male fertility in the adult (12). On the other hand, Dmrt1 binds to the Sox9 and Sox8 promoters and in Dmrt1(-/-) XY mutant gonads their expression was reduced (13). The interactions between Sox9, Sox8 and Dmrt1 shown in Fgf9 and on Wnt4. The paracrine function of Fgf9 is represented by a negative interaction towards Wnt4. For simplicity, we assumed that Foxl2 is capable of exerting its function without the support of its co-factor, the oestrogen receptor ERα. Finally, we considered that, unless the presence of Dmrt1 prevents it, Foxl2 is constitutively activated by 12.5 dpc in both sexes due to a developmentally programmed activator (see Section 2.3).
The positive interaction of Dmrt1 upon Sox9 was eliminated because it was found to be irrelevant for gonadal sex determination. Fig. 1 Wt1-Sf1 interaction. Wt1 is required to initiate Sf1 expression because, even at the earliest time when Sf1 expression can be observed at 9.5 dpc (14), Sf1 expression is absent in both XY and XX embryos mutant for Wt1 (15). Moreover, Wt1 binds to a 674 bp fragment in the Sf1 promoter causing its activation (15).

Wt1-Sry and Gata4-Sry interactions.
As previously mentioned, several factors have been identified that, when mutated, prevent the initiation of the male program because Sry expression is very much affected (reviewed in (16) and cites therein). Wt1 and Gata4 were chosen to represent these activators of Sry.

Sf1-Sry interaction.
In vitro transfection studies showed that Sf1 protein binds to and activates Sry promoter (17,18). Both XY and XX mice mutant for Sf1 develop female genitalia (19). Since the absence of either Sf1 or Wt1 causes male-to-female sex reversal, proper activation of Sry requires the presence of both Sf1 and Wt1.

Sry-Sox9 and Sf1-Sox9 interactions.
Sox9 is a direct target of Sry (20,21). Chromatin immune-precipitation assays revealed that Sry protein binds to the TESCO enhancer sequence in the Sox9 promoter (22). Up-regulation of Sox9 requires both Sry and Sf1 activators: Sf1 can activate TESCO in vitro although at low levels by itself, and Sry is able to activate TESCO in vitro only in the presence of Sf1 (22). Moreover, ubiquitous expression of Sry in XX gonads results in the initiation of Sox9 expression only in Sf1positive cells (23).

Sox9-Sox9 interaction.
Sox9 can bind the TESCO element of its own promoter via the same Sry binding sites (22). However, by itself, Sox9 has no activity in co-transfection assays, requiring the presence of Sf1 to activate TESCO (22). Indeed, a direct interaction between the N-terminal domain of Sox9 and the C-terminal domain of Sf1 has been described for the regulation of Amh promoter (24).

Sox9-Sf1 interaction.
In vitro transfection analysis indicated that Sox9 could activate Sf1 promoter in a reporter construct with increasing amounts of Sox9 (25).

Sox9-Sry interaction.
A negative feedback turns off Sry transcription, but only if the level of Sry protein is above a critical threshold required to initiate testis development (20). This is consistent with findings where Sry expression was prolonged in sex-reversed ovotestes in embryos carrying a weak allele of Sry (26). This negative feedback is mediated by Sox9 since the rapid decrease of Sry expression occurs just after Sox9 activation: Sry expression is not turned off at 12.5-13.5 dpc in XY gonads with low levels of Sox9 (10) or in XY gonads where Sox9 function is inactivated (27).
Fgf9-Wnt4 interaction. At 12.5 dpc, Wnt4 is up-regulated in Fgf9(-/-) XY gonads but not in wild type XY gonads. This suggests that Fgf9 is necessary for the down-regulation of Wnt4 in differentiating XY gonads (28). Moreover, treatment of XX gonads with exogenous Fgf9 suppresses Wnt4 normal expression (28).  Fig. 1 were based on this interpretation.

Logical model definition
Having defined the regulatory components and interactions to be included in the network, we are ready to define the logical model, i.e. specify the ranges taken by the variables associated to the components, their logical parameters, and finally the updating scheme to be considered for the model dynamics (see Section 1). This section justifies the requirement for multi-valued variables for some components, lists the logical parameters of our 1-cell network model and finally, discusses the temporal constraints found to ensure a final phenotype (i.e. in each condition (XX vs XY gonad), the reachability of the expected final state,).

Genes having more than one single functional level
Sf1 functional levels. In the bi-potential gonad of both sexes, the gene Sf1 is initially activated by Wt1 and, following Sry activation, becomes up-regulated in the XY but not in the XX gonad (25). We thus assumed that Sf1 initial expression corresponds to its first functional level, whereas its Sox9-dependent up-regulation corresponds to its second functional level. This second functional level is justified because Sox9 function requires Sf1 as co-factor and because Sox9 up-regulation produces more Sox9 protein, requiring a higher amount of Sf1.

Sox9 functional levels.
It was found that Sox9 auto-regulation needed to be operative already at its initial expression level in the bi-potential gonad for the model reproduce its sexual development. Hence, we assumed that Sox9 initial activation by Sf1 corresponds to its first functional level, and that subsequent expression of Sry brings Sox9 to its second functional level. This second level would be required to trigger the male development of the bi-potential gonad.

Fgf9 and Wnt4 functional levels. We assumed that the initial expression of both Fgf9
and Wnt4 (level 1) relates to their roles in cell proliferation during the formation and developmental plasticity of the bi-potential gonad. The specific up-regulation of Fgf9 in XY gonads after the increase of Sox9 expression (28) would constitute Fgf9 second functional level, involved in the inactivation of Wnt4-signalling pathway (28). Wnt4 second level would reflect its up-regulation in XX gonads, being involved in the inhibition of Sox9 auto-regulation (28).

Logical parameters
The values of the logical parameters define the behaviour of the regulatory network (see

Stable states analysis
GINsim enables the identification of all the stable states of a logical model. Using this feature, 27 stable patterns were identified (data not shown). However, this number greatly decreases if we restrict the input values to relevant combinations (see Supporting Table S3). These results suggest that the temporal signals are instrumental for the selection of the adopted differentiated state. Morever, as shown by the model simulations, the selection of a relevant initial condition leads to trajectories that recapitulate observed differentiation decisions by reaching the Sertoli or granulosa cell phenotypes.

Temporal constraints; defining priorities
Analysing the model dynamics ( Supplementary Fig. S2

XY gonad lacking Sry function developed into ovary (45), whereas gain-of-function
Sry mutation in XX gonad developed into testis, in agreement with experimental observations (46). (10,27,37) and XX gonad constitutively expressing Sox9 gave rise to testis, in agreement with experimental results (47). Moreover, model simulations predicted that Sox9 partial loss-of-function XY gonad would develop into ovary, whereas Sox9 partial gain-of-function XX gonad still gave rise to ovary.

Fgf9 loss-of-function caused XY gonad to develop into ovary, in agreement with
experimental observations (28,29,48,49 or intact (Fgf9 non-accessibility) gonads exposed to exogenous Fgf9. To this end, the initial state was that of the wild type XX gonad except that Fgf9=2, for the dissociated gonad (Fgf9=1 for the intact gonad); i.e., the high or low accessibility of Fgf9 was reflected by the high or low concentration of Fgf9 in the cell at the initial state. In the case of a dissociated gonad, the exogenous Fgf9 induced Sox9 expression, whereas such induction did not occur for an intact gonad (data not shown). These results provided a formal demonstration to the explanation of (28).

5.
Lack of Wnt4 function in XX gonad resulted in testis development (50), whereas Wnt4 ectopic expression in XY gonad resulted in ovary development, in agreement with experimental results (51). As previously mentioned, ß-catenin is the effector molecule of Wnt4 signalling pathway. Accordingly, mutations in ß-catenin are expected to mimic the mutations in Wnt4. This was indeed the case as XX gonad lacking ß-catenin function was predicted to develop as testis. Furthermore, XY gonad constitutively expressing ß-catenin developed as ovary, in agreement with experimental observations (52).

6.
It has been reported that XY gonads lacking Dmrt1 function initially develop as normal testes but become feminised later (36,43). The same effect has been reported when Dmrt1 function is removed from foetal or adult testes (13). The simulation of XY gonad lacking Dmrt1 function agreed with those experimental observations, in accordance with the role played by Dmrt1 in maintaining testis identity. In addition, gain-of-function Dmrt1 mutations were predicted to determine testis development of XX gonad. In this respect, it has been reported very recently that the expression of a conditional Dmrt1 transgene in the ovary silences the female sex-maintenance gene Foxl2 and consequently re-programmes juvenile and adult female granulosa cells into male Sertoli-like cells (53).

7.
It has been reported that XX gonads lacking Foxl2 function initially develop as normal ovaries but become masculinised later (33). Moreover, the loss of Foxl2 function in adult ovaries causes their trans-differentiation into testes (34). Simulating XX gonad lacking Foxl2 function reproduced those experimental observations, in accordance with the role of Foxl2 in maintaining ovarian. Moreover, gain-of-function Foxl2 mutations were predicted to determine the ovarian development of XY gonads.

Simulation of both XY and XX gonads carrying simultaneously loss-of-function Fgf9
and Wnt4 mutations develop as testes, in agreement with experimental results (30).

Mutations of the developmental temporal signals
The sexual phenotypes of the gonads carrying either loss-of-function (KO) or gain-offunction (GF) mutations in the genes encoding the developmental temporal signals AS, Gata4, IW and AF were simulated. The predicted phenotypes follow.  Our model accounts for these results (see complete network of Fig. S1). To this respect, the ensuing observations are relevant for explanation purposes. First Sox8 shows a significant degree of homology to Sox9 at the DNA-binding HMG box and its two flanking regions and a lower homology at the C-terminal region (56). Secondly, Sox8 is turn on in the male gonad by Sox9 after this is up-regulated in a Sry-dependent manner (57,58). Although Sox8 is a direct target of Sry (59), its main activator is Sox9 as the transcription of Sox8 directly depends on Sox9 expression (58). Thirdly, Sox8 activates the gene Amh, although this activation is weaker compared to Sox9-dependent activation (57), likely due to the differences between Sox8 and Sox9 in the C-terminal transactivation domain. And fourth, Sox9, throughout its C-terminal domain, interacts with ß-catenin resulting in their mutual destruction (32). These observations led us to assume that Sox8 and ß-catenin do interact also to form a complex that results in their mutual destruction, that is, that Sox8 and ß-catenin are engaged in a feedback loop.

Threshold number of Sertoli cells required for gonadal development into testis: the 2-cell network.
To address this question, two replicas of the 1-cell network of Fig. 1 were connected, defining a new network (Fig. 4A): the replica identified as network "c" represents the gonadal centre and the replica identified as network "p" represents the gonadal pole. In  Supporting Table S2).
In addition, to get the 2-cell network dynamics mimics the biological process, it was necessary to refine the temporal constraints previously defined, as a consequence of the inclusion of Fgf9-receptor (Fgf9r). Thus, faster transitions (priority class 1) were those increasing the levels of Wnt4, ß-catenin and Foxl2 and those decreasing the levels of Dmrt1, Fgf9r (in both cells), whereas all remaining transitions were considered slower.
The distinct temporal activation of Sry in central versus pole regions was modelled by assuming that AS, the developmental activator of Sry, exerts its function earlier in the central region -represented by AS_c-than in the pole region -represented by AS_p.
Hence, the initiation phase was subdivided into two sub-phases (t1 and t2), and the input signals were defined by: at t1, presence of Y, Wt1, Gata4, AS_c and absence of AS_p, AF and IW, and at t2, presence of Y, Wt1, Gata4, AS_c, AS_p and absence of AF and IW. In this scenario, the final state of t1 was taken as initial state of t2, whose final state was the initial state of the maintenance phase (t3).   Steroidogenic factor 1 Sry Sex-determining region Y Sox8 Sry-box 8 Sox9 Sry-box 9 Pgd2 Prostaglandin D2 synthase Fgf9 Fibroblast growth factor 9 Fgfr2 Fibroblast growth factor receptor 2 Wnt4 Wingless related MMTV integration site 4 Fzd G protein-coupled receptors, Class F frizzled Lrp Low density lipoprotein receptor-related protein Rspo1 R-sponding 1 Lgr4-6 Leucine-rich repeat-containing G-protein coupled receptor 4-6 b-cat Catenin beta-1 Dmrt1 Doublesex-and mab-3-related transcription factor 1 Foxl2 Forkhead-domain transcription factor L2 ERα oestrogen receptor alfa Amh Anti-Müllerian hormone Fst Follistatin