Assisted reproductive technologies (ARTs) represent a consolidated clinical practice, which have resulted in several million births since 1978 (Adamson et al., 2018) (International Federation of Fertility Societies’ Surveillance, 2019). However, the current ART protocols allow to use only a limited number of oocytes derived from antral follicles, while the large ovarian reserve, represented by the pool of early-stage follicles, remains a genetic patrimony that can be preserved by cryopreservation procedures but not managed (Martinez et al., 2017). In this context, the development of protocols aiming to obtain and fertilize mature oocytes from immature follicles grown outside the body could represent a useful strategy to recover the largest pool of ovarian gametes by promoting their in vitro growth and differentiation.

This follicle rescue approach could be applied as a possible future clinical strategy to preserve the ovarian reserve, as in the case of adult and prepubertal patients with cancer, where ovarian transplantation may expose to the risk of reintroducing malignant cells (Herta et al., 2018) (De Vos et al., 2014). Besides, it would benefit the veterinary medicine, providing a way to improve reproductive performance of species of zootechnical interest as well as to conserve the genetic inheritance of endangered animals.

Since the first attempt in the field of in vitro folliculogenesis (ivF) in 1996 (Eppig, 1996) (O’Brien et al., 2003), many efforts are made to set up new culture systems able to support in vitro growth of early-stage follicles toward competent oocytes (Laronda et al., 2017) (Xiao et al., 2015) (Xiao et al., 2017) (McLaughlin et al., 2018). The use of animal models to recapitulate ivF steps, driving to the production of fertilizable oocytes, has proven to be decisive for providing a knowledge basis and for developing validated methods with high translational potential for humans. Based on similarities in physiology and anatomy of the ovaries, folliculogenesis timing, and the follicle size (Barboni et al., 2011) (Bähr and Wolf, 2012) (Telfer and Zelinski, 2013), medium-sized mono-ovulatory mammals are commonly accepted as a translational model, and they are increasingly considered as being very relevant for human preimplantation reproductive research.

However, apart from the murine model where in vitro production of fertilizable oocytes has reached high levels of efficiency in terms of embryo development (Xu et al., 2006) (Hornick et al., 2013) (Hikabe et al., 2016), in medium-sized mammals, ivF remains still experimental. Indeed, a very low number of embryos produced from in vitro grown preantral follicles were reported in these models (de Figueiredo et al., 2018) even if several groups are working on bovine (Gupta et al., 2008) (Gupta and Nandi, 2012) (Antonino et al., 2019), porcine (Wu et al., 2001), caprine (Magalhães et al., 2011), ovine (Barboni et al., 2011) (Luz et al., 2012) (Arunakumari et al., 2010), and non-human primates (Xu et al., 2009) (Xu et al., 2011).

The difficulty to recapitulate in vitro the process of folliculogenesis in non-rodent animal models appeared to be related to the longer period required for follicle/oocyte growth, the greater dimension of antral follicles with competent oocytes, and the difficulty to mimic the environmentally favorable conditions to guarantee a synergic oocyte and somatic compartment development by preserving the tissue architecture (Rossetto et al., 2016).

Considering the large amount of data collected in vitro on the molecular mechanisms involved in the folliculogenesis among different species and the advances in in vitro follicle culture models, involving in vitro 2D and 3D culture approaches (Jones and Shikanov, 2019), the adoption of mathematical models might represent a valuable tool to organize the evidence collected to date by offering predictive models.

This study supports the use of a computational method based on network theory to identify the molecular events and the main factors sustaining ivF steps in mammals and to discover new molecular players in the ivF process to be targeted and/or exploited for therapeutic purposes, improving female fertility.

The ivf process is described by a scale-free, non-clustered, and non-hierarchical system network.

Up to half of the links within the network are referred to rodents (39%) and ruminants (ovine 17.2, caprine 3.9, and bovine 13.7%), whereas less than 9% comes from humans (Table 1).

Based on data collected in the WoS_MMivF database, it was possible to generate a biological network whose topological parameters were computed as reported in Supplementary File 3.

The network displayed a scale-free topology, according to the Barabási-Albert (BA) model. The statistical analysis of its topology (Table 2) demonstrated that it recognizes 641 nodes and 2089 links, with a very low clustering coefficient (0.076).

The main controllers of the network were identified: hub, bottleneck, and hub-bottleneck nodes.

The ivF network recognized as hyperconnected nodes (hubs) 5.9% of nodes (38 out of a total of 641) (Supplementary File 4), which were ranked based on node degree (Supplementary File 4). They include activating systemic hormones and ovarian activating factors (8 and 11 out of 38, respectively; overall 50% of hubs) more than transduction terminal events (4, 10, and 5%).

A Kernel density estimation (KDE) based on their clustering coefficient was then performed (see Supplementary File 4 for clustering coefficient value) in order to explore their distribution, showing that no node subpopulation is present (as reported in Supplementary File 5).

Furthermore, to assess the hub role in the ivF network stability control, a computational experiment was performed, carrying out subsequent cycles of hub removal (Supplementary File 6). The network topology was deeply affected by hub removal and collapsed upon two cycles of network attack. As comparison, the removal of the same number of nodes randomly identified [=random.between(min;max)] did not have detectable effects on network stability (data not shown).

In addition, seeking for controllers of the information flow within the network, 38 bottleneck nodes were identified (BN: Supplementary File 4) (Ordinelli et al., 2018) (Taraschi et al., 2020).

Finally, by intersecting the two subsets of nodes categorized as hubs and BN, 30 main nodes were identified and named hub.BN (Supplementary File 7).

The ivF network (Albert and Barabási, 2002) is characterized by a power-law distribution of node degree and by the absence of correlation between the node degree itself and the clustering coefficient in keeping with the Barabasi-Albert model, conferring to the network the following biologically relevant properties:

1) Robustness against random damage: when a random perturbation affects the network, it is very likely that only a scarcely linked node (i.e., a node belonging to the most frequent class of nodes) will be affected. Thus, the probability that a hub of the graph is affected remains very low, and in this model, it can be estimated to be about 5.9%.

2) Controllability: the small number of highly linked nodes implies that the whole system can be modulated with high efficiency by acting just on a few molecules, thus reducing the energetic cost and facilitating/accelerating the cell response.

3) Easy navigability: the virtual absence of clustering, together with the low values of the characteristic path length and of the average number of neighbors, confers to the network a typical structure of signaling networks. Interestingly, this low value of clustering coefficient also implies a little redundancy.

After having defined the topology, the computational analysis allowed to dissect the network to identify the ivF network flow (hub.BN). Indeed, the 2D KDE approach allowed to select those with a major modulatory role on the network information flow by identifying either cultural functional endpoints or in vitro follicle controllers such as systemic hormones and local factors (Matzuk et al., 2002).

More in detail, the priority of hub.BN was assigned by 2D KDE analysis to four events related to both follicle (primordial to primary follicle transition and preantral follicle growth) and oocyte development (oocyte meiotic competence). Furthermore, the analysis identified FSH as a key hub.BN; thus, reinforcing the idea of its central role as a controller of either the gonadotropin-dependent or independent phase of folliculogenesis (Holesh and Lord, 2018) (Casarini and Crépieux, 2019). Several pieces of scientific evidence collected mainly using in vitro studies, indeed, suggest the need to supplement FSH to also stimulate the early stage phases (from primary to later stages) (Hsueh et al., 2015) of folliculogenesis by stimulating either follicle growth or differentiation (Cecconi et al., 1999) (Barboni et al., 2011) (Cadoret et al., 2017) (Candelaria et al., 2020) (Cortvrindt et al., 1997) (Kreeger et al., 2005) (Wright et al., 1999). In addition, 7 hub.BN have been identified as main controllers of the network (activin A and LH, GDF-9, EGF, E2, Insulin, and IGF-1) (Orisaka et al., 2006) (Dong et al., 1996) (Edson et al., 2009).

The results of the network stratification analysis are in agreement with current knowledge on the physiological cross-talk between molecules regulating the inter- and intra-follicular communication. Indeed, in the input layer were identified hub.BN recapitulating the main modulatory factors, either hormones exerting a remote control on the ovary (LH, FSH, hCG, IGF1, and EGF) or molecules involved in intra-ovarian control of both somatic and germinal compartments (TGF-beta superfamily members and growth factors).

The processing layer is composed of hub.BN molecules belonging to PI3K/AKT, JAK/STAT, TGFbeta, and KL signaling pathways and the second messenger cAMP which transduce and amplify the actions of the previous endocrine and paracrine controllers (Li et al., 2008) (Fujihara et al., 2014) (John et al., 2009) (Imai et al., 2014).

Finally, the output layer recognized steroidogenic molecules (P4 and CYP19A1) which may be considered as endpoints of the pathways controlling ivF. The emerging picture defined by the computational stratification analysis thus supports the biological strength of the network recapitulating and organizing in three layers the physiological feedback orchestrating ovarian follicle development.

The biological strength of the network has been further confirmed by analyzing the reproductive phenotype of the available knock-out mice generated by using some ivF selected hub.BN. The analysis of the scientific evidence collected to date not only confirmed that the in vivo silencing of these controllers always promoted a negative impact on female fertility but also confirmed the interference specifically on pathways controlling the early phase of in vivo folliculogenesis, despite the different underlying mechanisms. STRING analysis was then performed to find new molecules and define their functional connections. Taking advantage of such an approach, molecules belonging to three effector categories have been identified as potential controllers of the early stage of the in vitro follicle development. Metabolic controllers involved in the signaling leading to follicle recruitment, apoptosis, and differentiation were identified. Indeed, STRING enriched the ivF network with mTOR and IGF1 signaling pathway molecules.

mTOR represents a molecular sensor for diverse environmental inputs including nutrients and growth factors and regulates various fundamental processes including cell growth, the metabolism, differentiation, and autophagy (Long et al., 2003). In folliculogenesis, some lines of evidence have shown that the IGF1-dependent activation (hub.BN role in the ivF network) of the PI3K/AKT signaling (hub role in the ivF network) and its downstream cognate mTOR complex sustains ovarian primordial follicle dormancy and activation, oocyte maintenance and activation, and GC proliferation and differentiation (Makker et al., 2014) (Guo and Yu, 2019) (Herta et al., 2018).

More in detail, mTORC1 plays a documented role in primordial follicle-oocyte bidirectional signaling. It was shown to activate the KIT receptor in oocytes through the KIT ligand (hub role in the ivF network), which triggers a PI3K/PTEN/AKT/forkhead box 3 (FOXO3) cascade and awakens the dormant oocytes (Zhao et al., 2018) (Liu et al., 2014). In the awakened oocyte, secretion of oocyte-specific growth factors such as BMP15 (hub.BN role in the ivF network) and GDF9 (hub.BN in the ivF network) further activates receptor serine kinases and downstream SMAD (mothers against DPP homolog 1 Drosophila) proteins including SMAD4 (BN role in the ivF network) in surrounding GC, leading to their growth and proliferation (Sanfins et al., 2018). A targeted deletion of the mTORC1 negative regulator TSC2 (STRING-enriched partner) in mouse oocytes results in prematurely follicular activation due to elevated mTORC1 activity in oocytes which in turn cause depletion of follicles in the early adulthood (Adhikari and Liu, 2009). Of note, a compensatory elevation of PI3K signaling was proposed to the reason for the unaffected follicular development observed in RPTOR (STRING-enriched partner) conditional KO in primordial and all subsequent oocyte stages (Gorre et al., 2014). Conversely, conditional KO mice for RPTOR in primordial follicle GC prevent the cell differentiation, and this arrests the dormant oocytes in their quiescent states, leading to the oocyte death age (Liu et al., 2014), indicating that the KIT/PI3K cascade in oocytes is indispensable for primordial follicle survival.

The key role of mTOR signaling was also supported by the computational identification of the downstream target of mTOR, such as FOXO1, FOXO3, and SIRT1 (STRING-enriched partners). Furthermore, transgenic mice for these molecules have provided evidence on their cross-talk in regulating the dynamics of the primordial follicle pool (Adhikari and Liu, 2009) (Adhikari et al., 2013) (Chen et al., 2015).

FOXO3 functions at the earliest stages of follicular growth as a suppressor of follicular overactivation, increasing the follicle reserves in the ovary in order to extend the reproductive period of females (Castrillon et al., 2003) (Shah et al., 2018) (Lee and Chang, 2019). Accordingly, FOXO3 has been found to be highly expressed in the nuclei of oocytes of primordial follicles, and its expression is downregulated in oocytes of primary and later-growing follicles, indicating that its downregulation in oocytes could be a prerequisite for the initiation of oocyte growth during follicular activation (Liu et al., 2007).

Furthermore, FOXO3−/− female mice exhibit a distinctive ovarian phenotype of global follicular activation, leading to oocyte death, early depletion of functional ovarian follicles, and secondary infertility (Castrillon et al., 2003) (Hosaka et al., 2004) (Lin et al., 2004).

Of note, selective depletion of FOXO1 and FOXO3 in mouse GC leads to an infertile phenotype characterized by metabolic changes and the production of factors that exerts potent negative feedback to prevent gene expression of pituitary FSH (hub.BN role in the ivF network). Decreased levels of serum FSH further restrict follicle growth and development, ultimately preventing ovulation. Besides, FOXO1/3 depletion alters the expression of genes involved in follicle growth and apoptosis, disrupting cell regulatory signals associated with the granulosa cell metabolism and follicle growth. These results reveal a novel ovarian-pituitary endocrine feedback loop preventing uncontrolled proliferation and/or premature differentiation of GC in follicles where apoptosis is impaired (Liu Z. et al., 2013).

Accordingly, the STRING-enriched molecule FOXO1 was shown to act as a silent guardian of follicle development. Indeed, it is a critical factor in promoting GC apoptosis to counteract the stimulatory FSH role, one of the selected ivF network hub.BN (Liu et al., 2009) (Fan et al., 2010) (Shen et al., 2014). More in detail, FSH was found to promote the expression of GC genes required for proliferation, survival, and estrogen synthesis by decreasing FOXO1, which negatively regulates proliferation and steroidogenesis (Liu et al., 2009) (Fan et al., 2010) (Shen et al., 2014) (Rosairo et al., 2008).

Moreover, the identification of SIRT1 as a potential partner of STRING enhances the role of the mTOR pathway in the control of in vitro early mammalian folliculogenesis. Indeed, it has been proposed as the molecular mediator of the calorie restriction-dependent prevention of follicular activation, leading to ovarian reserve preservation (Bordone et al., 2007) (Long et al., 2019). At a molecular level, the effect of SIRT1 on the ovary occurs through the deacetylation dependent-activation of FOXO3 and the suppression of the mTOR signaling (Long et al., 2019). The latter could be mediated by TSC2-SIRT1 interactions (Ghosh et al., 2010). These results are in good agreement with data showing that the SIRT1 activator (SRT1720) improves the follicle reserve and prolongs the ovarian lifespan of diet-induced obesity in female mice via activating SIRT1 and suppressing mTOR signaling (Zhou et al., 2014). Conversely, a recent work shows how SIRT1 can sustain the activation of mouse primordial follicles independent of its deacetylase activity. Specifically, SIRT1 was shown to trigger primordial follicle awakening by activating the PI3K/AKT-AKT and TSC1/2-mTOR signaling pathways. Indeed, its pharmacologically (resveratrol) and/or genetic-induced activation in mouse ovary cultures leads to SIRT1 to work as a transcription cofactor which modulates the expression levels of AKT, mTOR, and genes related to classic primordial follicle activation (Zhang et al., 2019). Interestingly, SIRT1 overactivation was shown to inhibit FOXO3 activity by promoting its exclusion form the nucleus through the AKT-dependent phosphorylation (Zhang et al., 2019), suggesting that SIRT1 activators might be used to efficiently activate the primordial follicles to be applied in in vitro activation protocols and/or to avoid uncontrolled follicular atresia characterizing POF.

The evidence of mTOR and SIRT1 pathways on follicular fate has been achieved by taking advantage of their in vivo and in vitro pharmacological modulation.

In vivo, the SIRT1 activation has been achieved by caloric restriction or the application of specific sirtuin activators or mTOR blockers. Caloric restriction was shown to induce the accumulation of SIRT1 in murine ovaries, which is associated with inhibition of primordial follicle activation and impairment of ovarian follicle development (Tatone et al., 2015) (Liu et al., 2015) (Long et al., 2019) mechanisms probably involved in mice delay of puberty, prolongation of reproductive lifespan, and prevention of age-associated infertility (Liu M. et al., 2013).

Various effects on folliculogenesis have been observed in vitro by using the SIRT1 activator, mTOR inhibitor, and resveratrol.

More in detail, it promoted growth in human ovarian follicles (Hao et al., 2018) and in ovine primordial follicles (Bezerra et al., 2018). The resveratrol action is addressed to regulate GC cell proliferation and survival (Han et al., 2017) (Ortega et al., 2012) (Schube et al., 2014) (Wong et al., 2010) as well as steroidogenesis (Qasem, 2020) (Morita et al., 2012) even if with contradictory results.

In late folliculogenesis, mTOR signaling is dependent on gonadotropin FSH modulation to manage processes such as oocyte maturation (Guo et al., 2018), ovarian somatic cell proliferation, and steroidogenesis (Palaniappan and Menon, 2012) and to intensify EGFR/RAS signaling (Fan et al., 2009) (Gloaguen et al., 2011) (Fan et al., 2012) (Wayne et al., 2007).

Interestingly, the STRING network enrichment identified the isoform KRAS as a potential interactor of the AKT1 and EGF hub.BN molecules with a high confidence score, suggesting its key role as a metabolic effector of gonadotropin stimulation. Its contribution in follicular development is corroborated using conditional knock-in mouse models in which the GC expresses a constitutively active form of KRAS (KrasG12D). These models show alterations in GC differentiation, proliferation, and apoptosis at early stage of folliculogenesis, thus impairing cell responses to gonadotropins and leading to premature ovarian failure (Fan et al., 2008) (Fan and Richards, 2010) (Bulun et al., 2019). This phenotype is similar to that of the mutant mouse model with EGFR signaling defects (Hsieh et al., 2007). The altered response of KRASG12D-expressing GC to gonadotropins appears to be related to low levels of FSHR and the inability of FSH to induce expression of LHCGR mRNA and therefore to the loss of the crucial LH-MAPK3/1 signaling pathways. This conclusion is supported by the reduced expression of specific genes known to be essential for COC expansion and ovulation (Richards, 2005). Moreover, mutations causing KRAS overactivity have been reported as markers of epithelial and endometrioid ovarian cancer (Ramalingam, 2016) (Dinulescu et al., 2005).

Although the effect of these energy sensors in folliculogenesis remains to be investigated in depth, some data have been reported in the literature.

The third category of molecules enriched from the STRING network was represented by several angiogenetic factors. Indeed, it must be considered that as the follicle growth, the mTOR-dependent metabolic signaling needs to be implemented with components of the vascular system that guarantees the correct trophic supply and spreading of precursors to properly complete folliculogenesis (Martelli et al., 2017).

The angiogenesis process assumes a key role in the contest of the reproductive function (Robinson et al., 2009) (Wulff et al., 2002). Accordingly, angiogenesis inhibition leads to the reduction of follicular growth, ovulation disruption, and drastic effects on the corpus luteum activity and development (Robinson et al., 2009). The large amount of factors is able to control angiogenesis and their spatio-temporal regulatory activities, suggesting that more than one factor might be useful for angiogenesis well-functioning associated with ovulation.

Specifically, functional interactions from STRING were unveiled for VEGF members (VEGFA, VEGFB, VEGFC, and VEGFD), VEGF receptors (NRP1, KDR, and FLT), and angiogenic factors (FGF1, FGF9, PGF, and FGFR).

Several pieces of evidence confirmed that both somatic and germinal compartments contribute to modulate the expression of VEGFA by modulating the blood vessel network during preantral follicle development. More in detail, small and middle pig preantral follicles seem to behave as autonomous recruited units, where the growth of the somatic compartment is always accompanied by the simultaneous activation of the endothelial cells (Martelli et al., 2017).

Differently, once the preantral follicles reach the late stage enclosing an almost fully growth oocyte, it starts to express high and stable levels of VEGFA, essential to maintain follicle angiogenesis in a steady status of activation (Martelli et al., 2017). Furthermore, high levels of VEGFA also characterize the preantral follicular structures at the stage of antrum formation. In this context, VEGFA has been supposed to increase microvessel permeability by stimulating plasma extravasations, thus allowing the accumulation of fluids within the differentiating follicular cavity (Isobe et al., 2005).

The strict correlation between somatic and vascular parameters may represent an indirect biological validation of the key role of VEGF members identified by ivF network STRING analysis depicting a synergic action between follicle compartments and blood vessel system components as a key event driving preantral follicle activation first and then sustaining the process of transition from preantral to early antral follicles (Martelli et al., 2017) (Danforth et al., 2003) (Abramovich et al., 2006).

These molecules exert their function through their receptors (Stouffer et al., 2001) and are tightly regulated by paracrine factors such as nitric oxide. Nitric oxide has been reported to mediate positive effects on follicle development and selection related to angiogenic events and play a modulatory role in the ovarian steroidogenesis (Basini and Grasselli, 2015).

Nevertheless, proangiogenetic factors such as FGF1, FGF9, and PGF have also been described as they can potentiate the effect of VEGF signaling by increasing the VEGF expression in theca cells of cattle (Gray et al., 1987) (Gospodarowicz et al., 1987) (Berisha et al., 2004).

In conclusion, the computational analysis of the ivF network has enabled us to select among several molecules adopted to date those that are the main spatio-temporal controllers of the early stage of follicle development. Taking advantage of the biological robustness of such a network, new molecules slightly deepened or unexplored for ivF purposes were identified by taking advantage of the STRING approach. Altogether, this evidence suggests the ivF network as a sounded system biology tool to be exploited for research, technological, and innovation aims with the final goal of designing new diagnostic and therapeutic strategies for female fertility.