Graphical abstract

   graphic file with name fx1.jpg
   [33]Open in a new tab

Highlights

     * •
       Differentiation of uniparental human ESCs into granulosa-like cells
       (GLCs)
     * •
       Androgenetic ESCs differentiate more efficiently to GLCs than
       parthenogenetic ESCs
     * •
       GLC differentiation of androgenetic ESCs depends on the imprinted
       gene IGF2
     * •
       Supplementation of IGF2 partly rescues GLC differentiation of
       parthenogenetic ESCs
     __________________________________________________________________

   In this report, Keshet and colleagues explore the genomic imprinting
   effect on ovarian development by differentiating biparental,
   parthenogenetic, and androgenetic human pluripotent stem cells into
   granulosa-like cells. They show that parthenogenetic cells display a
   reduced differentiation capacity, due to lack of expression of the
   paternal gene IGF2. These findings unravel an important contribution of
   the paternal genome to ovarian development.

Introduction

   Genomic imprinting is an epigenetic phenomenon in mammals, leading to a
   parent-of-origin dependent monoallelic expression of a group of genes
   that are required for proper embryonic development ([34]Barlow and
   Bartolomei, 2014; [35]Reik and Walter, 2001; [36]Tucci et al., 2019).
   Various genetic disorders arise from aberrations in parental
   imprinting, and disruption in the expression of imprinted genes plays
   important roles in many pathologies ([37]Peters, 2014).

   Dissecting the unique roles of the paternal vs. the maternal genome
   during development is extremely challenging due to the difficulty in
   uncoupling them. To this end, the use of uniparental embryos or
   embryonic stem cells (ESCs), which either have a strictly maternal
   (parthenogenetic or gynogenetic) or paternal (androgenetic) genome has
   contributed greatly to imprinting research.

   More than 3 decades ago, pivotal experiments in mice showed that
   uniparental embryos are not viable, thus implicating imprinting as the
   mammalian barrier against asexual reproduction ([38]McGrath and Solter,
   1984; [39]Surani and Barton, 1983; [40]Surani et al., 1984).
   Interestingly, these experiments also showed that while gynogenetic
   embryos had compromised extra-embryonic lineages, androgenetic
   conceptuses showed a compromised development of the embryo proper
   ([41]McGrath and Solter, 1984; [42]Surani and Barton, 1983; [43]Surani
   et al., 1984). Subsequent experiments with chimeric embryos composed of
   biparental and uniparental cells demonstrated different tissue
   distribution between parthenogenetic and androgenetic cells ([44]Barton
   et al., 1991; [45]Fundele et al., 1989, [46]1990; [47]Mann et al.,
   1990; [48]Nagy et al., 1989; [49]Thomson and Solter, 1988). Together,
   these observations suggested that there is a non-equivalent
   contribution from the two parental genomes to the development of
   various tissues.

   In humans, it is not possible to study embryonic development of
   chimeric or uniparental embryos in vivo due to obvious ethical reasons.
   Nonetheless, modeling tissue development through differentiation of
   uniparental human pluripotent stem cells (hPSCs) to different cell
   types can be utilized as a powerful alternative ([50]Stelzer et al.,
   2011). We have recently derived several novel human androgenetic and
   parthenogenetic ESC (aESCs/pESCs) lines. A global comparison of the
   methylome and transcriptome between these two cell types enabled the
   identification of new imprinted genes both within and outside of known
   imprinted regions ([51]Sagi et al., 2019).

   The importance of maintaining proper expression of imprinted genes
   during embryogenesis throughout adulthood is by now highly appreciated.
   Imprinted genes are known to regulate both behavioral and physiological
   processes including growth and survival, metabolism, the circadian
   clock, and neuronal function ([52]Peters, 2014). Nevertheless, while
   these various functionalities occur in the soma, the involvement of
   imprinted genes during the formation of the germline, especially in
   humans, remains elusive. This question is important from both
   evolutionary and reproductive therapy points of view. To address this
   issue, it is important to study the development of the germ cells as
   well as the somatic gonadal cells, which guide their differentiation.
   Indeed, if balanced imprinting is required for gonadal development, it
   adds a deeper layer of understanding to the mammalian barrier in
   asexual reproduction.

   Granulosa cells serve as somatic supporting cells to oocytes in the
   ovary. During mammalian fetal development, differentiation of granulosa
   cells from the gonadal primordium establishes the first essential step
   of female sex determination ([53]Rotgers et al., 2018). Indeed, these
   cells are not only vital for ovarian development, but are essential for
   ovarian function in adults ([54]Richards and Pangas, 2010). Thus, we
   reasoned that granulosa cells can be used as an attractive model for
   studying the development of the gonad. To explore whether genes that
   are expressed exclusively from one of the two parental genomes are
   essential for ovarian development, we compared the ability of
   biparental ESCs (biESCs), aESCs, and pESCs to differentiate into
   granulosa-like cells (GLCs). Remarkably, we observe an enhanced
   differentiation potential of aESCs compared with pESCs into GLC fate.
   Further analysis led us to identify IGF2 as the top-most upregulated
   paternally expressed gene (PEG) during the differentiation. Indeed,
   while knocking out IGF2 in aESCs results in their reduced
   differentiation capacity into GLCs, supplementation of recombinant IGF2
   to pESCs during the differentiation leads to a partial rescue in their
   differentiation efficiency. We thereby reveal an essential role for the
   paternal genome in the formation of a female gonad, which is at least
   partly explained by the PEG IGF2.

Results

Differentiation of biparental ESCs into GLCs

   To study the contribution of imprinted genes to the development of the
   human somatic reproductive tissue, we took advantage of a
   differentiation protocol, which enables the derivation of human GLCs
   from human ESCs ([55]Lan et al., 2013). This 12-day protocol starts
   with the formation of embryoid bodies (EBs), which are first guided
   into mesodermal fate and finally to GLCs, positive for the marker genes
   AMHR2, FSHR, CYP19A1, AMH, and FOXL2. Using this protocol, we
   differentiated biESCs into GLCs (biGLCs) ([56]Figure 1A). Around day 10
   of differentiation, large, vesicle-containing cells migrated out from
   the attached EBs ([57]Figure 1B). On day 12 of differentiation, we
   sorted AMHR2^+ biGLCs using fluorescence-activated cell sorting (FACS)
   for RNA extraction and performed RNA sequencing (RNA-seq). In agreement
   with previous reports ([58]Lan et al., 2013), day-12 biGLCs showed a
   consistent upregulation of AMHR2, FSHR, AMH, and FOXL2 granulosa
   markers and a downregulation of pluripotency markers compared with
   undifferentiated biESCs ([59]Figures 1C and [60]S1A). We noticed that
   sorting biGLCs led to very low yields and high degradation levels of
   RNA, possibly due to increased sensitivity of the differentiated cells
   to the sorting process. We therefore examined whether we could capture
   the transcriptional changes that occur throughout the differentiation
   in unsorted biGLCs samples. In agreement with the sorted samples,
   RNA-seq analysis of unsorted biGLCs showed the same upregulation and
   downregulation trend of different granulosa markers and pluripotency
   markers, respectively, albeit in a milder manner ([61]Figures 1D and
   [62]S1B). To further characterize the differentiation of biGLCs, we
   performed a differential expression analysis (DEA), which showed a
   significant upregulation of 1,467 genes and a downregulation of 554
   genes in unsorted biGLCs compared with undifferentiated cells
   ([63]Figure S1C and [64]Table S1). Indeed, when performing a gene set
   enrichment analysis (GSEA), one of the most enriched gene ontology (GO)
   terms was "urogenital system development," together with some
   mesodermal lineages, as expected from the GLC differentiation protocol
   ([65]Figure S1D). In further agreement, the "ovarian steroidogenesis"
   pathway was significantly upregulated when performing GSEA for Kyoto
   encyclopedia of genes and genomes (KEGG) categories, further supporting
   the efficiency of the differentiation ([66]Figure 1E). Recently, a
   detailed single-cell road map of human gonadal development during the
   first and second trimester was published ([67]Garcia-Alonso et al.,
   2022). The large number of cells analyzed in this study enabled the
   identification of different populations of granulosa cells along the
   granulosa cell differentiation trajectory. These populations included
   the first and second waves of pre-granulosa cells (preGC-I/preGC-II),
   together with follicular granulosa cells, which emerge from early/mid-8
   postconceptional week (PCW), and 17 PCW, respectively. Using this
   published dataset, we searched for differentially expressed genes that
   separate between the granulosa lineage (both pre- and follicular
   granulosa cells) and other lineages in the dataset (germ cells and
   other somatic lineages). This comparison yielded 240 genes that we coin
   "fetal granulosa markers" ([68]Table S2). Indeed, when performing a
   GSEA analysis on our biGLCs, upregulation of fetal granulosa markers is
   significantly enriched during the differentiation from biESCs
   ([69]Figure S1E), strengthening the ability of GLC to model granulosa
   cells during development. Next, we also performed a DEA to distinguish
   between different populations of fetal granulosa cells. This analysis
   resulted in 96, 70, and 268 markers that separate between preGC-I,
   preGC-II, and follicular granulosa cells, respectively ([70]Table S2).
   To try to position the GLCs along the granulosa cell developmental
   timeline, we compared the average expression of the top 10 markers for
   each fetal granulosa cell population in biESCs vs. biGLCs
   ([71]Figures 1F and [72]S1F). Interestingly, while markers for every
   fetal granulosa cell type were upregulated in GLCs compared with the
   undifferentiated state, the upregulation for preGC-I markers was the
   strongest, suggesting that GLCs correspond mainly to the first wave of
   pre-granulosa cells that emerge at 8 PCW ([73]Figure 1F). This result
   aligns with the relatively short differentiation protocol used in this
   study. Nevertheless, it should be noted that this similarity is
   determined only by a transcriptome analysis, and additional assays will
   be needed to measure the exact similarity of GLCs to in vivo preGLC-I.
   Together, these results suggest that biESCs can be efficiently
   differentiated into biGLCs, and that unsorted samples still capture the
   main transcriptional changes that occur during the differentiation
   process. We therefore continued the rest of the analysis using unsorted
   GLC samples.

Figure 1.

   [74]Figure 1
   [75]Open in a new tab

   Differentiation of biparental hESCs into GLCs leads to upregulation of
   granulosa specific markers

   (A) General scheme of GLCs differentiation protocol.

   (B) Microscopic images of day-10 GLCs. Scale bars, 50 μm.

   (C) Expression levels of different granulosa markers of sorted
   biparental GLCs (independent duplicates) vs. undifferentiated cells
   (independent triplicates [[76]Sagi et al., 2019], see [77]experimental
   procedures under “[78]external RNA-seq samples”). Bars represent mean ±
   SEM.

   (D) Expression levels of different granulosa markers of unsorted
   biparental GLCs (independent triplicates) vs. undifferentiated cells
   (independent duplicates). Bars represent mean ± SEM.

   (E) GSEA plots for GO categories that are enriched during
   differentiation to biparental GLCs and are likely relevant for
   granulosa cells.

   (F) Average expression (logCPM) levels of the 10 top-markers for
   different fetal granulosa cell types in unsorted biGLC vs. biESC cells.
   Bars represent mean ± SEM.

Biased differentiation of uniparental ESCs into GLCs

   To check whether differentiation of ESCs into GLCs is dependent on
   genes that are expressed exclusively from the paternal or maternal
   genome, we differentiated parthenogenetic or androgenetic ESCs into
   GLCs (pGLCs/aGLCs, respectively, three different cell lines each). RNA
   was isolated from each cell type and the transcriptome of biGLCs was
   compared to that of pGLCs and aGLCs ([79]Figure 2A and [80]experimental
   procedures). First, we focused on the 1,467 genes that were
   significantly upregulated during the differentiation of ESCs into
   biGLCs (referred to as "granulosa-associated genes"). The change in
   gene expression in the differentiation into biGLCs was then compared
   with the changes in gene expression from pESCs into pGLCs and from
   aESCs into aGLCs. Strikingly, while aGLCs upregulated
   granulosa-associated genes in similar levels to that of biGLCs, pGLCs
   failed to do so for most genes ([81]Figure 2B). We then performed a
   GO/KEGG pathways overrepresentation analysis for granulosa-associated
   genes, which were upregulated (logFC >1) during aGLCs differentiation
   compared with pGLCs differentiation (see [82]experimental procedures
   and [83]Table S3). The five top GO terms and KEGG pathways for these
   genes included "urogenital system development," "ovarian
   steroidogenesis," and "cholesterol metabolism," with "urogenital system
   development" being the top overrepresented GO term ([84]Figure 2C).
   Indeed, "reproductive system development" term was also significantly
   overrepresented, despite not being part of the top five categories
   ([85]Figure 2C). We also compared the expression changes of granulosa
   markers that showed consistent upregulation during differentiation of
   sorted and unsorted biGLCs (i.e., AMHR2 and FSHR) between biGLCs and
   aGLCs/pGLCs differentiation. Interestingly, this comparison showed that
   while aGLCs upregulated the granulosa marker AMHR2 at a higher level
   than biGLCs, pGLCs upregulated AMHR2 to a lower extent ([86]Figure 2D).
   Indeed, an initial FACS experiment performed with an antibody against
   AMHR2 showed a stronger signal in aGLCs compared with pGLCs
   ([87]Figure S2A). This staining was not optimal, however, enabling the
   identification of a subtle qualitative trend, rather than a
   quantitative result. This trend was also true for FSHR, albeit pGLCs
   had a negligible disadvantage in upregulating FSHR compared with
   biGLCs, suggesting that FSHR upregulation might be enhanced by some
   paternal genes but is not dependent on them ([88]Figure 2D). Together,
   these results suggest that the paternal genome has a greater
   contribution compared with the maternal one for GLC differentiation.

Figure 2.

   [89]Figure 2
   [90]Open in a new tab

   Androgenetic hESCs differentiate more efficiently to GLCs than
   parthenogenetic hESCs

   (A) Schematic diagram of the different cell types used for GLCs
   differentiation.

   (B) Comparison of androgenetic or parthenogenetic differentiation with
   biparental differentiation (aESC: independent duplicates; aGLCs:
   independent triplicates; pESCs: independent duplicates; pGLCs:
   independent triplicates), focusing only on genes that are upregulated
   during biparental differentiation into GLCs. p = 2.0e-27 by
   Mann-Whitney U test.

   (C) Overrepresentation analysis for genes that are upregulated during
   biparental differentiation and are also upregulated during androgenetic
   differentiation compared with parthenogenetic differentiation (with
   logFC >1).

   (D) Comparison between expression changes of granulosa cell markers
   during androgenetic/parthenogenetic differentiation compared with
   biparental differentiation.

IGF2 is the most upregulated imprinted gene during GLC differentiation

   Considering the advantage of aESCs to differentiate toward GLC fate, we
   analyzed the expression dynamics of imprinted genes during biGLC
   differentiation in search for candidates that might stand behind the
   bias. Interestingly, this analysis showed that the most upregulated
   imprinted gene during GLC differentiation is the PEG IGF2
   ([91]Figure 3A). In support of this observation, re-examination of the
   developing human ovary small-cell RNA-seq (scRNA-seq) from
   [92]Garcia-Alonso et al. (2022) showed that IGF2 is widely expressed in
   the developing ovarian somatic tissue ([93]Figure S2B). Specifically,
   in the granulosa lineage, IGF2 highest expression is observed in
   preGC-I (the cells to which GLCs seem to resemble the most). During the
   second wave of pre-granulosa cells, IGF2 levels seem to diminish and
   are elevated in follicular granulosa cells ([94]Figure S2B). This
   observation opens the possibility that IGF2 might serve as a paracrine
   as well as an autocrine factor during granulosa cell development. As
   expected, a comparison between the expression changes of imprinted
   genes that occur upon aGLC vs. pGLC differentiation showed that IGF2 is
   much more upregulated during aGLC differentiation ([95]Figures 3B and
   [96]S2C). To further confirm the imprinting status of IGF2 in our
   uniparental cell lines, we also analyzed an Infinium
   HumanMethylation450 BeadChips dataset ([97]Sagi et al., 2019) and
   calculated the H19/IGF2 imprinting control region (ICR) methylation
   levels in aESC and pESC ESC lines. As reference, we also included three
   sperm samples and biparental ESC lines ([98]Figures S2D and S2E and
   [99]experimental procedures). While the sperm samples presented
   hyper-methylation and the biparental cell lines showed ∼50% methylation
   levels, our aESC and pESC lines showed hyper- and hypo-methylation
   levels, respectively, along the H19 ICR ([100]Figures S2D and S2E). In
   addition, we examined the methylation status of KvDMR1, a maternally
   methylated differentially methylated region (that should therefore be
   hypomethylated on the paternal allele), which like the H19 ICR, sits in
   the Beckwith-Wiedemann syndrome imprinted region on chromosome 11. As
   expected, KvDMR1 presented a mirror phenotype compared with the H19
   ICR, with aESC and pESC lines showing hypo- and hyper-methylation
   levels, respectively ([101]Figures S2D and S2E). These results suggest
   that the uniparental cell lines maintain their parent-of-origin
   imprinting signature in the H19/IGF2 imprinted region. We therefore
   chose to focus on the role of the imprinted gene IGF2 during
   differentiation of ESCs into GLCs. Insulin-like growth factor 2 (IGF2)
   is a member of the insulin/insulin-like growth factors (IGFs) system,
   which also includes insulin and insulin-like growth factor 1 (IGF1).
   IGFs and insulin exert their physiological effects by binding and
   activating two tyrosine kinase receptors, namely the type-1
   insulin-like growth factor receptor (IGF1R) and the insulin receptor
   (INSR) ([102]Gallagher and LeRoith, 2010). While IGF2 can activate both
   the INSR and IGF1R, insulin can only activate the INSR ([103]Neirijnck
   et al., 2019). We thus hypothesized that IGF2 may mediate IGF1R
   signaling during GLCs differentiation. In agreement with this
   hypothesis, IGF2 showed much higher expression levels compared with
   IGF1 in biGLCs ([104]Figure 3C). Moreover, both INSR and IGF1R were
   highly expressed in biGLCs, with IGF1R showing the higher expression
   between the two ([105]Figure 3C). This observation makes sense
   considering that IGF2, and not IGF1, is known to be the main IGF1R
   signaling mediator in the human ovary ([106]El-Roeiy et al., 1993). In
   agreement with IGF2 being the only imprinted IGF family member, neither
   IGF1, IGF1R, nor INSR showed any notable differences in expression
   among androgenetic, parthenogenetic, or biparental GLCs
   ([107]Figure S2F). These results, combined with the fact that the
   differentiation media, which is supplemented with serum replacement,
   contains insulin but not IGFs, suggest that IGF2 acts as the main IGF
   in GLCs and could potentially mediate IGF1R signaling in these cells.

Figure 3.

   [108]Figure 3
   [109]Open in a new tab

   Different imprinted genes are upregulated during GLCs differentiation

   (A) Volcano plot representation for imprinted (colored) and unimprinted
   (gray) gene expression changes during biparental GLCs differentiation.

   (B) Volcano plot representation for imprinted (colored) and unimprinted
   (gray) gene expression changes during androgenetic vs. parthenogenetic
   differentiation.

   (C) Expression levels of IGF2/IGF1 and INSR/IGF1R in biparental
   granulosa cells. Bars represent the mean ± SEM.

The paternally imprinted IGF2 is essential for efficient GLC differentiation

   To explore whether IGF2 expression contributes to the higher propensity
   of aESCs to differentiate toward GLCs, we differentiated an aESC line
   harboring a knockout (KO) mutation for IGF2 (termed
   aESC/aGLCs^IGF2_KO), which we previously generated using CRISPR-Cas9 in
   combination with a single-guide RNA (sgRNA) targeting a constitutively
   expressed exon ([110]Sagi et al., 2019). Strikingly, RNA-seq analysis
   of aGLCs^IGF2_KO showed that they failed to upregulate most of the
   granulosa-associated genes compared with aGLCs ([111]Figure 4A).
   Interestingly, GSEA for GO terms and KEGG pathways that were associated
   with biGLCs and were also enriched during differentiation of aGLCs
   compared with pGLCs, showed that these categories were significantly
   depleted in aGLCs^IGF2_KO compared with aGLCs ([112]Figure 4B). Indeed,
   when comparing the full transcriptome of aGLCs^IGF2_KO with GLCs from
   different parental backgrounds, KO of IGF2 has a profound effect on the
   differentiation into GLCs ([113]Figure S3A). These results suggest that
   IGF2 strongly contributes to the advantage of aESCs, compared with
   pESCs, in differentiating toward GLCs. We also analyzed the expression
   levels of AMHR2, which was upregulated in higher levels during aGLCs
   and downregulated during pGLCs differentiation compared with biGLCs
   differentiation. This analysis showed a downregulation of the granulosa
   marker AMHR2 in aGLCs^IGF2_KO compared with aGLCs ([114]Figure 4C).
   Together, these results demonstrate an overall reduction in the
   efficiency of GLCs differentiation in the absence of IGF2, suggesting
   an essential role for IGF2 during the development of granulosa cells.
   Because IGF2 works in a paracrine/autocrine manner, a fascinating
   question is whether pGLCs differentiation could be enhanced if
   performed in the presence of exogenous IGF2. To explore this
   possibility, we supplemented pESCs with recombinant IGF2 throughout the
   differentiation protocol (termed pGLCs^+IGF2). Indeed, RNA-seq analysis
   of pGLCs^+IGF2 showed a significant enrichment of the most depleted
   categories in aGLCs^IGF2_KO, i.e., "reproductive system development"
   and "urogenital system development" categories, and an upregulation of
   AMHR2 in pGLCs^+IGF2 compared with pGLCs, essentially mirroring the
   results observed in aGLCs^IGF2_KO ([115]Figures 4D and 4E). To
   conclude, we observed an increased contribution of the paternal genome
   toward differentiation into GLCs, which could partly be linked to the
   PEG IGF2.

Figure 4.

   [116]Figure 4
   [117]Open in a new tab

   IGF2 is essential for GLC differentiation

   (A) Expression levels of genes upregulated during biGLCs
   differentiation in aESCs, aGLCs, or aGLCs^IGF2_KO (aGLCs^IGF2_KO:
   independent duplicates).

   (B) Relevant categories from GSEA based on aGLCs^IGF2_KO vs. WT aGLCs
   differential expression analysis. Dashed line represents significance
   threshold (p < 0.05).

   (C) Expression levels of AMHR2 for aGLCs^IGF2_KO vs. WT aGLCs. Bars
   represent mean ± SEM.

   (D) Relevant categories from GSEA based on pGLCs^+IGF2 (independent
   duplicates) vs. pGLCs (independent duplicates) differential expression
   analysis. Dashed line represents significance threshold (p < 0.05).

   (E) Expression levels of AMHR2 for pGLCs^+IGF2 vs. pGLCs. Bars
   represent mean ± SEM. NES, normalized enrichment score.

   Next, we searched for candidate TFs, which potentially mediate IGF2
   contribution during GLCs differentiation. To this end, we performed an
   overrepresentation analysis of TF binding sites for the promotors of
   granulosa-associated genes that are more upregulated during aGLCs
   compared with pGLCs differentiation and are also upregulated (logFC >1)
   in aGLCs compared with aGLCs^IGF2_KO ([118]Table S4). This analysis
   identified six TFs, namely IGLV5-37, TCF7, GREB1, HES2, ZNF596, and
   ZNF92, that potentially regulate these genes ([119]Figures S3B and
   S3C). Interestingly, one of the most enriched TF binding sites was for
   GREB1, which is a chromatin-bound estrogen receptor (ER) coactivator
   that is essential for ER-mediated transcription ([120]Mohammed et al.,
   2013). Indeed, examination of the expression levels for the six TFs in
   different tissues as obtained from the GTEx database showed that GREB1
   is highly enriched in the ovary compared with other tissues
   ([121]Figure S3D). This observation is even strengthened, considering
   that when looking on the scRNA-seq dataset from [122]Garcia-Alonso
   et al. (2022), GREB1 expression seems to be almost exclusively
   expressed in pre-granulosa and granulosa cells ([123]Figure S3E). These
   results possibly link between IGF2 signaling and GREB1-mediated
   transcription activation during granulosa cells development.

Discussion

   Understanding the unique contribution of each parental genome to the
   development and maintenance of the various tissues in the human body is
   of great importance for basic research as well as for the better
   understanding of the genetic disorders that are associated with
   aberrations in imprinting ([124]Peters, 2014). Imprinted genes have
   been extensively studied in various physiological and behavioral
   contexts. Nonetheless, whether the requirement of sexual reproduction
   in mammals, which is mediated by imprinted genes, extends to the
   formation of the sex organs remained unknown.

   In this work we have differentiated ESC lines from different parental
   backgrounds into GLCs and compared their differentiation capacity
   through their transcriptome. We show that upon differentiation, biGLCs
   upregulate granulosa markers and activate genetic programs that are
   associated with reproductive system development and ovarian steroid
   production. Strikingly, while aESCs differentiate to GLCs with similar
   efficiencies compared with biESCs ([125]Figure 2B), pESCs fail to
   upregulate most of the granulosa-associated genes to the same levels as
   biGLCs. Moreover, we show that the granulosa-associated genes that are
   upregulated more in aGLCs than pGLCs differentiation, are highly
   enriched in terms that are specific to ovarian development and
   steroidogenesis.

   This finding is especially interesting considering that previous
   studies have shown that parthenogenetic cells, but not androgenetic
   cells, are able to contribute to the germline in chimeric mice
   ([126]Fundele et al., 1990; [127]Nagy et al., 1989). More recently,
   these observations were strengthened by a report of a successful
   differentiation of parthenogenetic ESCs into functional oocytes, which
   upon fertilization led to the birth of live and fertile mice ([128]Tian
   et al., 2021). Importantly, this does not mean that germ cell
   development in mice is not dependent on PEGs, as gonadal supporting
   somatic biparental cells in the chimeric mice and in the in vitro
   differentiation protocol (which depends on culturing the primordial
   germ cells with somatic ovarian cells taken from a wild-type (WT) fetal
   ovary) could still rely on PEGs.

   Imprinting is conserved between mice and humans, but most imprinted
   genes are not shared between them ([129]Tucci et al., 2019), and the
   function of some of the shared genes is not identical (as discussed
   below regarding IGF2). Thus, whether there is also a barrier for
   androgenetic cells to differentiate toward primordial and functional
   germ cells in humans, remains to be investigated. Nonetheless, if such
   a barrier exists, it may indicate that while a maternal genome is
   needed for oocyte differentiation, the paternal genome is needed for
   the development of the supporting somatic cells in the gonad, thus
   creating a requirement for both parental genomes in female gonadal
   development.

   In search of PEGs that could potentially contribute to the
   differentiation bias of the paternal genome, we identified IGF2 as the
   most upregulated PEG during GLCs differentiation. Importantly, we were
   able to show that aGLCs^IGF2_KO fail to upregulate most of the
   granulosa-associated genes and reproductive system development genes
   compared with WT aGLCs. Furthermore, addition of recombinant IGF2
   throughout the differentiation of pESCs remarkably led to an
   upregulation of reproductive system development genes and the granulosa
   marker AMHR2 in pGLCs. The findings regarding IGF2 are intriguing
   considering that the insulin/IGF system is known to play important
   roles in gonadal function ([130]Neirijnck et al., 2019). Specifically,
   mice carrying a targeted disruption of the igf1r in granulosa cells are
   sterile, present small ovaries, and are devoid of antral follicles
   ([131]Baumgarten et al., 2017). Of note, granulosa IGF1R signaling in
   mice is likely mediated mainly by IGF1 and not IGF2, given that while
   igf1-null female mutants are infertile, paternally transmitted igf2
   disruption leads to the birth of small but fertile mice ([132]Baker,
   1996; [133]DeChiara et al., 1990). Moreover, IGF1 and not IGF2 is the
   primary IGF in the murine ovary ([134]Hernandez et al., 1989). In the
   human ovary, however, IGF2 is the primary IGF, and is most highly
   expressed by granulosa cells ([135]El-Roeiy et al., 1993). In the adult
   human ovary, IGF2 was shown to mediate the steroidogenic and
   growth-promoting actions of follicle-stimulating hormone (FSH) on
   pre-antral follicles during the menstrual cycle, through IGF1R-Akt
   signaling ([136]Baumgarten et al., 2014; [137]Yuan and Giudice, 1999).
   Specifically, FSH was shown to upregulate IGF2, which then synergizes
   with FSH to promote the proliferation and steroid production of
   granulosa cells in an Akt-dependent manner ([138]Baumgarten et al.,
   2015). However, while IGF2 has been strongly associated with the growth
   and differentiation of human adult follicles, whether it is essential
   for the initial development of granulosa cells, as modeled by our
   system, was previously unknown. Our findings thus add another layer to
   the known functions of IGF2 during human follicle maturation and
   highlight the importance of IGF2 already in the initial formation of
   granulosa cells. Interestingly, some patients with Beckwith-Wiedemann
   syndrome, which show biallelic expression of IGF2, were reported to
   suffer from ovarian hyperplasia, demonstrating the dosage sensitivity
   for IGF2 expression in the human ovary ([139]Mahat et al., 2019). Taken
   together, our results highlight the divergent roles of conserved
   imprinted genes such as IGF2 between mouse and human, considering that
   IGF1R signaling in the mouse ovary is primarily mediated by the
   non-imprinted gene igf1 and not igf2 ([140]Baker, 1996; [141]DeChiara
   et al., 1990). Moreover, we believe our findings should be taken into
   consideration when differentiating ESCs into GLCs, as IGF2 loss of
   imprinting is a known phenomenon that can happen when culturing ESCs
   ([142]Bar et al., 2017; [143]Keshet and Benvenisty, 2021), thus leading
   to a creation of a model that does not fully recapitulate its in vivo
   counterpart.

   Nevertheless, our study has some limitations. First, while the
   differentiation was strong enough to be captured in the transcriptome,
   it is important to keep in mind that because we did not sort the GLCs,
   the transcriptome is expected to represent some noise coming from a
   certain degree of differentiation heterogeneity in the samples.
   Moreover, despite the differences we observe between aGLCs and pGLCs,
   it is yet to be determined how these transcriptomic differences
   translate to the functionality of granulosa cells.

   Finally, the construction of reconstituted ovaries from somatic cells
   and germ cells that will all be differentiated from pESCs isolated from
   patient oocytes hold great promise for reproductive medicine. Our
   findings thus may be adopted in such future protocols to mitigate the
   lack of PEGs during GLCs differentiation, through addition of
   recombinant IGF2.

Experimental procedures

Resource availability

Corresponding authors

   Further information and requests for resources and reagents should be
   directed to and will be fulfilled by the corresponding authors, Nissim
   Benvenisty ([144]nissimb@mail.huji.ac.il), Gal Keshet
   ([145]gal.cleitman@mail.huji.ac.il), and Talia-Eldar-Geva
   ([146]gevat@szmc.org.il)

Materials availability

   This study did not generate new unique reagents.

Cell lines

   Human ESC lines used in this study included three previously derived
   androgenetic cell lines (aES1, aES3, aES5) ([147]Sagi et al., 2019),
   three previously derived parthenogenetic cell lines (SWAPs4, pES6,
   pES10) ([148]Paull et al., 2013; [149]Sagi et al., 2016), and three
   previously derived biparental female lines (WA09, CSES7, CSES10)
   ([150]Thomson et al., 1998; [151]Lavon et al., 2008). For culture
   conditions and experiment repetitions see [152]supplemental
   experimental procedures.

Differentiation of ESCs into GLCs

   Differentiation into GLCs fate was performed as previously described
   ([153]Lan et al., 2013) with slight modifications. See
   [154]supplemental experimental procedures.

Flow cytometry

   For sorting of GLCs, differentiated GLCs were dissociated using TrypLE
   Select (Thermo Fisher Scientific) into single cells. The cells were
   then quenched with FACS buffer (10% FBS in PBS), centrifuged, and
   resuspended in FACS buffer containing mouse anti-human AMHR2 (1 μg per
   10^6 cells, Abcam ab64762) antibody for 30 min on ice. Next, the cells
   were washed twice with FACS buffer followed by 30-min incubation with
   an Alexa Flour 488-conjugated Donkey anti-Mouse antibody (1:800, Abcam
   ab150109). The cells were then washed twice with FACS buffer. For each
   sort, a small number of GLCs were set aside and were only stained with
   the secondary antibody to serve as negative controls. At the end of
   each staining, cells were filtered through a 70-μm cell strainer and
   analyzed by flow cytometry (BD Biosciences FACSAria III) and FlowJo
   software (FlowJo LLC). For the flow cytometry experiment in
   [155]Figure S2A, differentiated GLCs were dissociated using TrypLE
   Select (Thermo Fisher Scientific) into single cells. The cells were
   then quenched with FACS buffer (10% FBS in PBS), centrifuged, and
   resuspended in FACS buffer containing rabbit anti-human AMHR2
   cy3-conjugated antibody (4 μg per 10^6 cells, antibodies-online,
   ABIN1702944) for 30 min on ice. Next, the cells were washed twice with
   FACS buffer, resuspended in FACS buffer, filtered through a 70-μm cell
   strainer and analyzed by flow cytometry (BD Biosciences FACSAria III)
   and FlowJo software (FlowJo LLC).

Generation and analysis of IGF2-knockout aESC lines

   The IGF2^KO aESC line (aES3) was previously obtained in the
   lab ([156]Sagi et al., 2019). For the IGF2^KO differentiation
   experiment, two independent replicates of aES3 IGF2^KO were used, and
   differentiation was performed using the same conditions as in other
   cell lines.

RNA isolation and RNA-seq

   For high-throughput RNA-seq experiments, differentiated and
   undifferentiated cells were harvested at the end of differentiation
   assays. Total RNA was isolated using RNeasy Mini Kit (QIAGEN) and the
   mRNA fraction of total RNA was enriched by pulldown of poly(A)-RNA. RNA
   sequencing libraries were generated using TruSeq RNA Library Prep Kit
   (Illumina) according to the manufacturer’s protocol and sequenced using
   Illumina NextSeq 500 with 75 bp single end reads. For transcriptome
   analysis, reads were mapped to the GRCh38 reference genome using STAR
   ([157]Dobin et al., 2013). Because the hESCs samples were grown on
   mouse feeder cells, all the samples (including the GLC samples) were
   subjected to the XenofilteR software ([158]Kluin et al., 2018) with the
   MM_threshold paremeter set to 15 in order to remove reads originating
   from mouse contamination. Count table was created using featureCounts
   ([159]Liao et al., 2014). The sorted biGLCs library was un-stranded.
   The undifferentiated controls for the unsorted GLCs and all unsorted
   GLCs were prepared with minus stranded library preparation kit.

External RNA-seq samples

   Because the sorted biGLCs library was un-stranded, three biESC samples
   that were prepared with an un-stranded library were downloaded from the
   SRA database to serve as controls for the sorted biGLCs ([160]Sagi
   et al., 2019) (accession number SRA: SRR7186863, SRR7186864,
   SRR7186853). The samples were of cell lines CU-ES4, CU-ES5, and NYSCF2.

Differential expression analysis

   In RNA-seq experiments, EdgeR package in R v3.32.0 ([161]Robinson
   et al., 2010) was used for DEA, using trimmed mean of M-values (TMM)
   method for normalization and the generalized linear model approach for
   model testing ([162]Dunn and Smyth, 2018; [163]Robinson and Oshlack,
   2010). For detailed description of the normalization, DEA model design
   and pathway enrichment analysis see [164]supplemental experimental
   procedures.

Single-cell RNA-seq of fetal granulosa cells analysis

   To extract general and specific fetal granulosa cell markers, the
   single-cell dataset named "Human female" was downloaded
   ([165]www.reproductivecellatlas.org) and subsampled to 100K cells using
   scanpy library in Python
   ([166]https://scanpy.readthedocs.io/en/stable/) for further analysis.
   The dataset was then converted into a Seurat object and the rest of the
   analysis was performed with Seurat package V4.2.0 in R
   ([167]https://satijalab.org/seurat/; see [168]supplemental experimental
   procedures).

Analysis of the methylation array

   Raw methylation data were downloaded from GEO under accession number
   GEO: [169]GSE114679 and analyzed using bioconductor’s minfi package
   (v.4.2). Genome browser style methylation representation was performed
   with bioconductor’s Gviz package (v.4.2). For detailed description of
   the analysis, see [170]supplemental experimental procedures.

Statistical analysis

   Hierarchical clustering was performed using the pheatmap
   ([171]https://cran.r-project.org/web/packages/pheatmap/index.html).
   Scatterplots, boxplots, and bar plots were generated in R using the
   ggplot2 package ([172]Wickham, 2016).

Author contributions

   G.K., S.B., N.B., and T.E.-G. conceived the study and designed the
   experiments. G.K. and O.Y. performed the experiments with a significant
   contribution from S.B. G.K. evaluated the data and performed the
   bioinformatic analysis. G.K. wrote the manuscript with input from S.B.,
   R.S.-G., N.B., and T.E.-G. N.B. and T.E.-G. supervised the study.

Acknowledgments