{'fraction_non_alphanumeric': 0.03315904723600125, 'fraction_numerical': 0.025427894713321714, 'mean_word_length': 5.101799672786766, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 3, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 1, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Single-cell roadmap of human gonadal development

Gonadal development is a complex process that involves sex determination followed by divergent maturation into either testes or ovaries1. Historically, limited tissue accessibility, a lack of reliable in vitro models and critical differences between humans and mice have hampered our knowledge of human gonadogenesis, despite its importance in gonadal conditions and infertility. Here, we generated a comprehensive map of first- and second-trimester human gonads using a combination of single-cell and spatial transcriptomics, chromatin accessibility assays and fluorescent microscopy. We extracted human-specific regulatory programmes that control the development of germline and somatic cell lineages by profiling equivalent developmental stages in mice. In both species, we define the somatic cell states present at the time of sex specification, including the bipotent early supporting population that, in males, upregulates the testis-determining factor SRY and sPAX8s, a gonadal lineage located at the gonadal–mesonephric interface. In females, we resolve the cellular and molecular events that give rise to the first and second waves of granulosa cells that compartmentalize the developing ovary to modulate germ cell differentiation. In males, we identify human SIGLEC15+ and TREM2+ fetal testicular macrophages, which signal to somatic cells outside and inside the developing testis cords, respectively. This study provides a comprehensive spatiotemporal map of human and mouse gonadal differentiation, which can guide in vitro gonadogenesis.


25
Summary: Gonadal development is a complex process that involves sex determination followed by divergent maturation into ovaries or testes. Historically, limited tissue accessibility and lack of reliable in vitro models have impeded our understanding of human gonadogenesis, despite its importance in gonadal pathologies and infertility. Here, we generated a comprehensive map of first-and second-trimester gonadal development using a combination 30 of single-cell and spatial transcriptomics, chromatin accessibility assays and imaging. Using this approach, we identified novel transcription factors and cell states in human germ and supporting cell lineages. We compared them with other mammalian species and found primatespecific regulatory programmes. Our data identified cell context-specific interactions shaping sex specification and development of human germ cells. We defined a novel bipotent 35 progenitor cell (LGR5+, TSPAN8+) in late embryos that can differentiate into early Sertoli in males or pre-granulosa cells in females. In fetal ovaries, we defined two subsets of pregranulosa cells supporting germ-cell differentiation and distributed across the corticomedullary axis. We also found a subset of developing granulosa cells appearing during the second trimester of pregnancy that is involved in follicular assembly. In fetal testes, we defined 40 a novel supporting population (sPAX8 cells) located at the poles of the developing testis cords. We also found two tissue-resident myeloid populations that we named microglia-like and SIGLEC15+ fetal testicular macrophages. This study provides an unprecedented spatiotemporal map of human gonadal differentiation that can be utilised as a blueprint for in vitro gametogenesis. 45 Gonadal development is a complex process requiring tight coordination of multiple cell lineages in time and space, some of which originate outside the gonads. Reproductive health of adults depends on these early developmental stages. Alterations in gonadal development 50 cause Differences in Sex Development (DSD), conditions with discrepancies between chromosomal, gonadal and phenotypic sex that are linked to infertility and, to a lesser extent, cancer 12 . In addition, gonadal dysfunction is associated with infertility, cancer and polycystic ovary syndrome, the most prevalent endocrinopathy among women of reproductive age 3, 4 . Limited access to fetal tissue, lack of in vitro models, and critical differences between humans 55 and other mammals have hampered our knowledge of human gonadogenesis 5 . Understanding this process will establish the foundation necessary for in vitro gametogenesis [6][7][8][9] , improve diagnosis of gonadal disorders 10 , and lead to design of effective therapies.
Primordial germ cells (PGCs) are the precursors of female and male gametes that transmit 60 genetic information to the next generation. Specification of PGCs occurs during gastrulation at ~2.5 weeks after fertilisation (post-conceptional weeks, PCW); thereafter PGCs migrate to the genital ridge 3-5 weeks after fertilisation 5, 11 . Sex determination is initiated in the somatic lineage by the action of SRY, the sex-determining gene located on the Y chromosome 12 . SRY expression leads to differentiation of Sertoli cells (supporting lineage), testosterone-producing 65 fetal Leydig cells (steroidogenic lineage) and other non-steroidogenic interstitial cells in XY gonads 13 . In XX gonads, somatic cells differentiate into pre-granulosa (supporting lineage) via Wnt4/Rspo1/β-catenin signaling and into ovarian stromal interstitial cells 14, 15 . Sex determination in the somatic compartment instructs germ cell differentiation as prespermatogonia in XY or oogonia and oocytes in XX gonads. Male germ cells enter mitotic 70 arrest, with germ cells and pre-Sertoli cells organising into cord-like structures surrounded by Leydig and other interstitial cells. At puberty, the testicular cords differentiate into seminiferous tubules where the male gametes reside. In contrast, female germ cells enter an asynchronous transition from mitosis to meiosis after ~10 PCW. Initially, oogonia develop in clusters named germ-cell nests that begin to break apart at ~17 PCW, when granulosa cells 75 surround primary oocytes to form primordial follicles that further differentiate into mature follicles at puberty.
Single-cell transcriptomics analyses on gonadal development have started to shed light on this process in humans [16][17][18][19] . These studies describe the main transcriptional changes in germ and 80 testicular somatic cell differentiation. However, holistic studies integrating multiple regulatory layers in space and time at the single-cell level are still missing. In addition, although much is known about development at puberty 20,21 , a detailed map of the anatomical organisation of the testis and ovary at a cellular level during fetal development is lacking. Understanding the close communication between germ and somatic cells that shapes sex specification and 85 differentiation is important, as well as identifying additional somatic cell lineages, such as tissue-resident macrophages, known to play a role in gonadogenesis in mice 22 .
Coupled analysis of single-cell and spatial transcriptomic atlases of tissues have laid the groundwork to understanding human development 23,24 , mapping disease 25, 26 and engineering 90 organoids 27, 28 . In addition, emerging studies of chromatin accessibility are revealing the first regulatory networks that have the potential to fine-map active transcription factors (TFs) during cell specification and associate genetic variants with cell type specificity [29][30][31] . Here, we profile 43 samples of human fetal gonadal and extragonadal tissue from the first and second trimester using a multi-omics approach that combines single-cell and spatial transcriptomics, single-cell 95 open chromatin profiling, and imaging. We compare our regulatory profiles to other mammals, gaining insight into the unique cell states and regulatory programmes activated in humans. We identify novel somatic cells: a bipotent supporting progenitor in the early gonad; two subsets of pre-granulosa and a subset of developing granulosa cells in the developing ovary; a supporting cell population plumbing the poles of the forming cords in the developing testis; 100 and two embryonic-derived testicular macrophages. We propose functions for the novel subsets according to their cellular neighbourhood and identify their potential implications in distinct gonadal conditions. 105 Results A single-cell multiomics and spatial map of developing gonads and extragonadal tissue We acquired gonad and extragonadal tissue from 43 embryos/fetuses in the first and second trimester of gestation (Carnegie Stage, CS16 -around 6 PCW-to 21 PCW) ( Fig. 1a-b). Samples 110 were classified as female (n=23) or male (n=20) based on karyotype and sex-specific gene expression (Supplementary Table 1). We used multiple genomics methods: i) single-cell RNA sequencing (scRNA-seq), ii) spatial transcriptomics (10X Genomics Visium), iii) singlecell accessible chromatin sequencing (scATAC-seq) and iv) combined single nuclei RNA and ATAC sequencing (snRNA-seq/scATAC-seq). Gonadal and extragonadal tissues were 115 processed in parallel and loaded in the same channel of the 10x Genomics chip.
Although we find typical endothelial, perivascular and erythroid cell populations in both gonadal and extragonadal tissue, mesenchymal cells differ based on location. We leveraged 135 Visium to spatially resolve transcripts in mesenchymal cells that reside inside and outside the gonadal tissue. We find that the expression of GATA4, LHX9 and ARX is exclusive to the gonads while expression of GATA2 and NR2F1 is restricted to extragonadal tissue ( Fig. 1f-g,  Supplementary Fig. 3a). Based on spatial information, and the expression of known and novel markers, we further annotated the mesenchymal clusters ( Supplementary Fig. 2b). We refined 140 our annotations by mapping the cells back into the tissue using cell2location 33 ( Supplementary Fig. 3b). Within the gonadal mesenchymal cells we find testis-specific interstitial cells in males, including fetal Leydig cells, ovary-specific stromal cells in females  and a population located at the dorsal mesentery shared by both sexes.  145  To gain insight into the chromatin landscape that shapes developing gonadal cells, we profiled  67,202 and 30,276 cells for female and male gonads, respectively, using scATAC-seq  (Supplementary Table 1-2). Batch effects were corrected using Harmony 32 and cell populations were annotated leveraging scRNA-seq data from matched individuals 34 (Fig. 1hi, Supplementary Fig. 4a-d). Cell identity assignment was validated by snRNA-seq data from 150 the combined snRNA-seq/snATAC-seq profiling (Supplementary Fig. 4e). We annotated a total of 17 populations in females and 15 in males and used chromatin accessibility peaks to: i) identify TF motifs that are differentially accessible in the distinct cell lineages 35 and ii) define genomic regulatory cis-co-accessibility networks (CCANs) specific to each cell type using Cicero 36 (Supplementary Fig. 4f-g, Supplementary Table 4). Our scATAC-seq data 155 confirms that accessibility of GATA4, LHX9, ARX, GATA2 and NR2F1 clearly differs between gonadal and extragonadal mesenchymal populations in both females and males ( Supplementary Fig. 4h-i).
Furthermore, we made use of our multiomics map to assign cell type specificity to the genes 160 that are frequently mutated in Differences in Sex Development (DSD) 1 ( Supplementary Fig.  3c-d). We find expression and motif activity of DSD-associated genes enriched in pregranulosa (EMX2, FOXL2), Sertoli (SOX9/10) and germ cells (DMRT1). Additional cell typespecific genes are described in further sections of our manuscript as we delve into each cell lineage. 165
PGCs are found at all stages of development, but their proportion decreases with fetal age, representing 100% of total germ cells at 6-7 PCW and dropping to 12% and 40% at 17-21 PCW in females 175 and males, respectively. In females, the STRA8 surge, which marks the initiation of meiosis 37 , appears at 11 PCW. Oocytes are apparent by 14 PCW and increase in frequency as they develop, accounting for 30% of female germ cells by 17-21 PCW. In males, pre-spermatogonia appear at 11 PCW and increase to nearly 60% by 17-21 PCW (Fig. 2b). 180 We used spatial transcriptomics and cell2location 33 to map cell states to the tissue architecture. The fetal ovarian cortex remains PGC rich at least until 17 PCW (latest stage analysed by Visium) (Fig. 2e). The ovarian medulla evolves to embrace a heterogeneous composition of states including mature germ cells and retinoic acid signaling-responsive oogonia at 11 PCW, meiotic oogonia at 14 PCW and oocytes at 17 PCW (Fig. 2e, Supplementary Fig. 5a). Female 185 germ cells differentiation across the forming cortico-medullary axis was validated by singlemolecule fluorescence in situ hybridisation (smFISH) with RNAscope probes (Fig. 2f). In contrast, PGCs in males are located in the forming testicular cords (Fig. 2f, Supplementary  Fig. 5b). 190 We then measured the regulatory dynamics of the TFs mediating germ-cell differentiation by taking into account their expression and activity, the latter measured by (i) the expression of consensus targets; and/or (ii) the chromatin accessibility of their binding sites in the germ cell scATAC-seq data (Fig. 2g, Supplementary Fig. 5c- Fig. 2g). Overall, there is a drastic change in the TFs and chromatin landscape when transitioning from PGCs to pre-spermatogonia. Conversely, the process of oogenesis is accompanied by more gradual changes in TF and chromatin accessibility dynamics (Supplementary Fig. 5f). 200 We find novel TFs mediating germ cell differentiation in humans and compare their expression in other mammalian species by re-analysing publicly available scRNA-seq datasets of developing testes (E10 to E16) and ovaries (E10 to E16 and E11.5 to P5) from mouse, as well as developing ovaries from macaque (E84 and E116) 38-40 ( Supplementary Fig. 5g-j, 205 Supplementary Table 6). Murine scRNA-seq data included pre-and postnatal samples to accommodate the different dynamics of germ cell differentiation and follicle development 41 .
One of the TFs, SOX4, is expressed in PGCs and spermatogonia but is downregulated during oogenesis (Fig. 2h). This expression pattern is conserved in mouse and macaque. The transition from PGCs to pre-spermatogonia involves the activation of the novel regulators EGR4, KLF7, 210 TWIST1 and NFIL3, displaying the same expression pattern in mice.
Fetal oocyte differentiation is more complex than its male counterpart, as it involves committing to the oogenic fate and transitioning from mitosis to meiosis. At the time of the STRA8 surge, we find activation of ZGLP1, the oogenic TF recently described in mice 42 . The 215 homeobox factors HOXA are upregulated in humans, a result that has not been previously described to our knowledge (Fig. 2g). The activation of these factors is not conserved in macaque and mice (Fig. 2h). At the meiotic stage, the TFs DMRTC2 and ZNF711, previously described in mice, are active. We also identify another conserved DMRT member, DMRTB1, which is upregulated in macaque and, to a lesser extent, in mouse oogonia (Fig. 2g-h). In 220 oocytes, we find activation of TP63, TSC22D1, conserved in macaque and mouse, and ZHX3, conserved in macaque but not mouse (Fig. 2g-h). TP63 is a gene associated with ovarian insufficiency 43 , showing the relevance of our TF roadmap to human germ cell development.
Human bipotent supporting cells are defined by co-expression of LGR5 and

TSPAN8
Gonadal supporting cells guide mammalian germ cell differentiation and maturation. In mice, supporting cells arise in, at least, two waves 39,44,45 . During the first wave, bipotent progenitor cells derived from the coelomic epithelium differentiate into pre-granulosa cells (females) or alternatively, if the Y-linked gene SRY is expressed, into Sertoli cells (males). This first wave 230 of supporting cells determines the sex specification of the gonad, which will in turn trigger the subsequent developmental waves. In humans, the bipotent cells giving rise to the first wave of male and female supporting cells have not yet been defined. Amongst the cells of the supporting lineage, we find a new population, sLGR5, which expresses known stem cells markers (eg. LGR5, TSPAN8) and other gonadal TFs (eg. GATA4) (Fig. 3a-c, Supplementary Table 7). 235 sLGR5 is enriched in the late stages of embryo development CS16-CS22 (~6-8 PCW) and over-expresses SRY in males (Fig. 3b-c; Supplementary Fig. 6a, Supplementary Table 7).
The expression profile of sLGR5 points at these cells as the bipotent progenitor that differentiate into Sertoli cells or the first wave of pre-granulosa cells. We then set out to validate its spatial location in developing testes and ovaries at early stages (CS19-CS20) using multiplexed RNA Scope ( Fig. 3e; Supplementary Fig. 6c). sLGR5 is mainly located in the 245 medulla, co-existing with the first wave of pre-granulosa cells (OSR1+) in females, and in close contact to developing Sertoli cells (SOX9+, LGR5-) in males. Therefore, the spatiotemporal dynamics of sLGR5 is in keeping with this population being the bipotent progenitor.
By 8 PCW, there are increased numbers of sex-specific supporting cells in males and females (Fig. 3b). In ovaries, three main subsets within the developing granulosa lineage (FOXL2+, 250 IRX3+) are present (see below). In testes, an early Sertoli population (Sertoli WFDC2) enriched between 8-10 PCW is connected to more advanced Sertoli cells (Fig. 3a-b). Sertoli cells have a unique extracellular matrix (ECM) composition facilitating interactions with germ cells to form the testicular cords (Fig. 3f, Supplementary Table 8). They express integrin alpha-6/beta-1 and NCAM1 that interact with junction adhesion molecules (JAM) JAM2/JAM3, as 255 well as other ECM proteins such as LAMC1, collagen IX (COL9A1/COL9A3) and AGRN. In the assembling testicular cords, specific signals triggered by the expression of SOX9 in the supporting cells orchestrate specification to the male lineage (Fig. 3f). Sertoli-to-germ cell signaling in the testicular cords is restricted to Inhibin β (INHA/INHBB) and TGFB1, known regulators of spermatogenesis. Indeed, cognate receptors for Inhibin β (BMPR2 and ACVR2A) 260 and TGFB1 (TGBR1-R2 complex) are upregulated in spermatogonia. WNT5B, which can activate the non-canonical WNT pathway in Sertoli cells, was one spermatogonia-specific ligand we identified. Sertoli cells also express oxytocin receptor, OXTR, which binds to its cognate ligand, OXT, expressed by all germ cells. In females, WNT signaling dominates in the distinct pre-granulosa subsets. In contrast, WNT ligands in the coelomic epithelium (WNT2B, 265 SOSTDC1, WNT5A) are repressed in male supporting cells, with the exception of the noncanonical WNT ligand WNT6 (Supplementary Fig. 6d). In line with this, Sertoli cells do not express AXIN2 , a well-defined target of the WNT canonical pathway 48 ( Supplementary Fig.  6d).
Although PAX8 expression is characteristic of Mullerian and Wolffian ductal epithelium 49 , both sPAX8 subsets express lower levels of EPCAM and do not cluster with ductal epithelial 280 cells when analysed in the same manifold ( Supplementary Fig. 7a-d). Therefore, sPAX8 is a new supporting population, distinct from the ductal epithelial PAX8+ cells.
We used smFISH to characterise sPAX8 in male and female gonads. In the early samples (CS19 and CS20), PAX8 expression is localised at the gonadal-mesonephric interface, between the 285 mesonephric tubules and the gonadal medulla (Fig.3g). Gonadal PAX8+ cells are EPCAM low , similar to the other gonadal supporting and germ cells. This is different from the EPCAM high expression pattern characteristic of the epithelial cells in the mesonephric tubules. The proximity of sPAX8 (EPCAM low , PAX8+) and mesonephric tubules (EPCAM high , PAX8+) in early stages, forming a continuum, suggests that epithelial-mesenchymal transition (EMT) 290 could be taking place as cells migrate into the gonad (Fig.3g). In line with this, sPAX8 expresses the TF SNAI2 involved in EMT 50 , as well as the detachment-induced protein PERP 51 , the suppressor of anoikis (apoptosis induced by detachment) NTRK2 52 , and the protease PLAU, all involved in controlling cell adhesion and apoptosis (Fig. 3c). We find differences between the location of sPAX8 in males and females in later stages (late 8 PCW onwards) (Fig.  295 3h, Supplementary Fig. 7e-f). In developing testes, sPAX8 cells are in the medulla, at the poles of the developing testicular cords, a location consistent with the rete testis ( Fig. 3h-i,  Supplementary Fig. 7e). In developing human ovaries, few sPAX8 cells are observed, maintaining their location at the ovarian-metanephric interface, near the hilum ( Supplementary Fig. 7f). Indeed, the presence of rudimentary rete ovarii has been described 300 at this interface in early development, but degenerates in later stages 53 .
The gonadal-mesonephric interface, where the sPAX8 is located, is a site of extensive tissue remodelling. At this interface, the movement of cells and their products within the gonads and formation of the rete testis are controlled (Fig. 3i). Accordingly, sPAX8 cells display a unique 305 transcriptional pattern of morphogens, suggesting a structural and supporting role for this population during development (Fig. 3h). Early sPAX8e expresses CXCL12/14 and its receptor CXCR4; the latter is also expressed by endothelial, epithelial and sLGR5, suggesting a chemotactic role for these populations (Fig. 3j-k). Male sPAX8m expresses NRP2, which binds VEGF proteins expressed by endothelial cells and SEMAB/C on the epithelial cells ( Fig. 3j-k).

310
Male sPAX8m also expresses additional molecules such as the somatostatin SST and IGFBP3, whose receptors are broadly expressed in Sertoli, pre-spermatogonia (TMEM219) and epithelial cells from the mesonephric duct (SSTR2) (Fig. 3j-k).
We then aligned our single-cell signatures with publicly available mouse female gonad scRNA-315 seq data using a support vector machine (SVM) model trained on human supporting cells ( Supplementary Fig. 7g-h, Supplementary Table 9). We find a probable murine sPAX8 cell population that shares expression of TBX1, TBX2 and LYPD1 with its human counterpart and also distinctly expresses other markers, such as Pdgfgl and Aldh1a3, which are absent in humans ( Supplementary Fig. 7i). In mice, sPAX8 are also located inside the gonads close to 320 the interface with extragonadal tissue (Supplementary Fig. 7j). These results suggest main features of sPAX8 cells are conserved in mammals.
TFs defining the two waves of medullary and cortical pre-granulosa cells in humans
Early preGC-I are located in the medulla, while preGC-IIa/b mapped to the outer cortex, subjacent to the ovarian surface epithelium (LRRN4+, KLK11+). preGC-IIc and developing granulosa are located in the transition zone between the outer cortex and the medulla. This region in the inner cortex will eventually become the location of the primordial follicles after 350 ~17 PCW. The spatial relationship among female supporting populations is reflected in their global chromatin accessibility patterns ( Fig. 4c-d, Supplementary Fig. 8b-d, Supplementary Table 7). All pre-granulosa subsets share a similar co-accessibility profile as the ovarian surface epithelium, with the exception of early medullary pre-granulosa (preGC-I). Thus, as previously described in mice, there appear to be two waves of pre-granulosa cells (medullary 355 and cortical) during human ovarian development 39 .
We then aligned the female supporting populations with publicly available female gonadal scRNAseq data from macaque and mouse using a support vector machine (SVM) model (Supplementary Fig. 8e-h, Supplementary Table 9). Overall, humans and macaques share 360 similar pre-granulosa cell populations, but, intriguingly, we could not establish direct correspondence with mouse pre-granulosa cells. We also could not establish a match of our bipotent supporting cell sLGR5 in mice.
LGR5 expression is restricted to bipotent progenitors in humans, while in mice LGR5 is expressed by the second wave of pre-granulosa cells, which are derived from the ovarian surface ( Supplementary Fig. 8i). This suggests that, despite the 365 presence of similar developmental waves of pre-granulosa cells in mice and humans, the regulatory networks giving rise to them are distinct.
We therefore explored the unique TF modules regulating human pre-granulosa development integrating transcriptomics and chromatin accessibility, as we did with the germ cells (Fig. 4e, 370 Supplementary Table 10). Some of the TFs reflect distinct lineage relationships of preGCs. The first wave of medullary preGC-I activates the master regulator of the granulosa lineage, FOXL2 54 , and the novel TFs STAT1, OSR1, ELF3 and MAFB, which are shared with the sLGR5 progenitor in females. The second wave of pre-GCs (preGC-II) expresses the previously described ARX, LHX2/9 55 and the novel EMX2 and ASCL1, shared with the ovarian 375 surface epithelium. Other TFs reflect environmental and functional differences. For example, preGC-I and preGC-IIc express WNT-induced TFs (FOXO1, FOXP1, HIF1 or FOXJ3), suggesting there is a higher WNT environment deeper in the ovary. In support of this, the WNT modulators SFRP1/2 are highly expressed in the ovarian surface and absent in any other cell state in the medulla ( Supplementary Fig. 6d). The developing granulosa cells express TFs 380 previously associated with germline nest breakdown and primordial follicle formation in mice, including NOTCH-induced TFs (HEY2 and HES4) and PEG3, as well as novel TFs such as TSHZ1, SOX5 or TEF ( Fig. 4a and 4e).
In humans, the numbers and quality of primordial follicles formed during fetal life determines 385 the reproductive lifespan and fertility of women in adulthood 56 . Polycystic ovary syndrome (PCOS) is the commonest cause of anovulation in women of reproductive age and is associated with disordered growth and development of the follicles 57 . In order to explore the link between genetics and fetal origin of adult PCOS, we aligned cell type-specific co-accessible regions from intragonadal female cells with polymorphisms associated with PCOS risk from the 390 GWAS Catalog 58 . Out of 43 PCOS-associated variants exceeding genome-wide significance, 30 map to our co-accessible regions, with some regions containing multiple variants. The majority of these regulatory regions are enriched in preGC-II, developing granulosa and ovarian surface epithelium, but to a lesser extent in preGC-I, indicating that the cellular target of PCOS could be supporting cells in the cortical and transition zones (Fig. 4f). In addition, 395 fine mapping of one variant, rs8043701-A, reveals that it might affect the co-accessibility pattern of TOX3, which is upregulated in preGC-II and developing granulosa cells ( Fig. 4g-h).

Crosstalk between germ and supporting cells shapes ovarian architecture
Our spatial analysis reveals a gradient of pre-granulosa and germ cells in the forming ovary. 400 To systematically map gonadal cell-cell communication in the ovarian microenvironments we updated CellPhoneDB by doubling the number of manually curated interactions. We also included interactions mediated by small molecules such as steroid hormones. To do this, we considered expression of the last enzyme of the metabolic pathway synthesizing the small molecule (Supplementary Table 11). In addition, we developed CellSign, a novel module of 405 CellPhoneDB v4, that links receptors with their downstream TFs, as assessed by manual curation. TFs are used as sensors of pathway activities, allowing us to identify sets of cell-cell interactions that affect downstream signaling predicted in the previous sections with a high degree of certainty ( Fig. 5a-b). 410 In the outer cortex, primordial germ cells co-exist with cortical preGC (preGC-IIa and b), coelomic epithelium and the ovarian surface. CellPhoneDB reveals that cortical supporting cells are a major source of paracrine factors for the PGCs, such as BMP4 (coelEpi) and KITLG (preGC-II), both required for the chemotaxis and maintenance of PGCs (Fig. 5c). STAT1 and STAT3, downstream of KIT, are accordingly active in PGCs (Fig. 5d). We also detect novel 415 interactions between supporting cortical cells and PGCs mediated by ephrins, EFNB2/EFNB3 (coelomic epithelium), with their downstream effector ELK1 activated in PGCs (Fig. 5c-d).
PreGC-I in the medulla as well as preGC-IIc and developing granulosa cells in the transition zone express molecules that mediate oocyte differentiation (Fig. 5f). preGC-I upregulate the enzymes involved in the production of estrogen, HSD17B1 and CYP19A1, which could bind 425 the estrogen receptor ESR1 expressed by oogonia STRA8+ (Fig. 5c). preGC-IIc expresses the enzyme ALDH1A2 involved in the synthesis of retinoic acid, whose product will bind to the retinoic acid receptor RARA and transporter CRABP2 expressed by oogonia STRA8+ (Fig. 5c).
Around 17 PCW, granulosa cells in the transition zone surround germ cells and begin to form 435 follicles, a process that involves bidirectional communication between the two cell types. For example, developing granulosa cells express NOTCH2/3 receptors, which are inhibited via DLL3 by PGCs, regulated via DLK1 by oogonia, and finally activated via JAG1 by oocytes to induce folliculogenesis in females (Fig. 5c). Other factors expressed by oocytes are GDF9, key for folliculogenesis, and enzymes involved in histamine and steroid hormone biosynthesis. 440 Their cognate receptors are exclusively expressed by developing granulosa cells. Finally, we find a unique expression pattern of extracellular matrix proteins in the oocytes and developing granulosa that will likely participate in the formation of the follicle (Fig. 5c). For example, oocytes in the medulla upregulate ITGA4, ITGA11 and E-cadherin (CDH1) to attach to the developing granulosa cells via collagen IV (COL4A1-4) and integrin ɑ-2/β-1 (ITGA2/ITGB1) 445 (Fig. 5e).
Finally, we investigated if the main roles of the pre-granulosa lineage subsets are preserved in mice, despite the transcriptomic differences we have shown between rodents and primates. We plotted the main identified ligand/receptor interactions in the medullary and cortical pre-450 granulosa subsets from the mouse dataset. Despite major differences between murine and human pre-granulosa cells, the first-wave of murine medullary pre-granulosa cells expresses enzymes involved in estrogen synthesis (eg. HSD17B1) while the second wave of murine cortical pre-granulosa cells expresses the oogenic factor BMP2 ( Supplementary Fig. 8j). This suggests that the location and main functional modules that define the first and second wave of 455 pre-granulosa cells in primates and rodents are conserved.

Two subsets of tissue-resident macrophages in the developing testis
Tissue-resident immune cells are thought to be important in gonadal development and function 460 22,61 . To further characterize the immune compartment of the fetal gonads, we sorted cells from 11 samples using the pan-leukocyte marker CD45 ( Supplementary Fig. 9a Fig. 6a; Supplementary Fig. 9b-d,  465 Supplementary Table 12). Myeloid cells are the most abundant during the first and second trimester, accounting for over 50% of CD45+ cells in the gonads (Supplementary Fig. 9e).
As expected, most fetal macrophages share the characteristic tissue-repair signature (Fig. 6h).
In the absence of scRNA-seq data from fetal bone, SIGLEC15+ macrophages are exclusively 490 present in developing testes. Microglia-like testicular macrophages cluster together with a similar subset in skin, suggesting these unusual immune cells might have a broader role beyond their function in the brain.
Using smFISH, we confirmed the presence of tissue-repair macrophages in developing gonads 495 of both sexes and further investigated the anatomical location of microglia-like and SIGLEC15+ macrophages in the fetal testes. Tissue-repair macrophages in the developing testes and ovaries are interstitial, co-existing with other interstitial and mesenchymal lineages ( Supplementary Fig. 10a). SIGLEC15+ macrophages are found outside the testicular cords, in close proximity to endothelial cells ( Fig. 6i; Supplementary Fig. 10b-c). Microglia-like 500 macrophages are usually found inside the testicular cords, with the potential to contact Sertoli and germ cells ( Fig. 6j; Supplementary Fig. 10d-e).
Together with the expression of immunomodulatory molecules (IL10, TGFB1, ENTPD1), a 515 potential function for microglia-like macrophages is in supporting an immunoregulatory environment inside the developing cords to avoid inflammation or oxidative stress that could damage maturing germ cells (Fig. 6l).

Discussion
Generating a cell atlas of gonadal development is essential to define the cell states and pathways involved in sex specification and divergent differentiation of ovaries versus testes. Our multiomics map, profiling more than one third of a million cells, will improve 525 understanding of gonadal conditions including DSD, infertility and cancer. Studying gonadal development is challenging, because it is a highly dynamic process that involves tight coordination of multiple cell lineages, including PGCs and the formation of oogonia and spermatogonia in females and males, respectively. Here, we have used single-cell transcriptomics, chromatin accessibility and spatial transcriptomics methods to generate a 530 comprehensive dataset of gonadal and extragonadal cells throughout embryonic and fetal development browsable online at reproductivecellatlas.org (user: rca, password: $Uj6mPXA). We discover novel supporting and immune cells, compare them to other mammalian species and thereby reveal primate-specific regulatory programmes. We identify the key interactions between somatic and germ cells that orchestrate the formation of male and female gonadal 535 structures during fetal development, and use this information to describe the potential function of new cell types.
Building on previous work profiling the transcriptome of human germ-cells 5,78,79 , our joint analysis of transcriptome and chromatin accessibility reveals novel TFs involved in germ cell 540 differentiation. Our multiomics analysis allowed us to prioritise TFs that are both expressed and active. Our multi-species comparison highlights a conserved TF network in primates with multiple differences from their murine counterparts, in agreement with previous murine-human comparisons of primordial germ cells 5 and spermatogonia 80 . Our work also reinforces the relevance of development in postnatal disease. For example, our analysis highlights the 545 expression of TP63 in oocytes, a gene associated with ovarian insufficiency 43 .
When developing germ cells enter the gonadal ridge (~5 PCW), sex-specific programmes are activated in somatic cells of the developing testes and ovaries (~6 PCW). We define for the first time in humans a bipotent supporting progenitor expressing unique markers (TSPAN8+, 550 LGR5hi, GATA4+) between 6-8 PCW. SRY expression peaks in this population in males. This bipotent population is likely to be derived from the coelomic epithelium 81 and gives rise to either the first wave of pre-granulosa cells or, if SRY is expressed, to Sertoli cells. A bipotent supporting progenitor has been previously described in mice 12 , although the murine markers (WT1+, NR5A1+) are broadly expressed by all supporting populations in humans. 555 The formation of testicular cords is a major event, indicating the start of male sex specification. We have identified essential interactions between Sertoli and germ cells that orchestrate differentiation and sex specification. Formation of the cords critically depends on testicular vascular and mesenchymal cell migration from the adjacent mesonephros to surround Sertoli-560 germ cell aggregates 77,82,83 . This process has to occur in a well-organised manner. We identify a supporting sPAX8 population (WNT6+, GATA4+) that is located in the gonadal-mesonephric border, peaks at 8 PCW, and is conserved in mice. After 13 PCW, this population is virtually absent in females and remains close to the developing cords in males. sPAX8 has a unique expression of chemokines and paracrine factors indicating this population has a role in 565 structural organisation by creating gradients of morphogens. Specifically, sPAX8 is located in one of the poles of the forming cords, suggesting it is involved in the formation of the rete testis, a complex network of tubules, into which the testis cords anastomose 84,85 . Indeed, PAX8+ cells are present in the adult rete testis 49 , the site of initial seminiferous tubule's fluid reabsorption. Altogether, this population may have a structural role organising the endothelial, 570 mesenchymal and epithelial cells that will form and connect the seminiferous tubules, rete testis, efferent ductules and the head of the epididymis.
In addition, we find that the developing testes harbor two unusual tissue-resident macrophages that are likely to participate in the testicular cord formation, and that we named microglia-like 575 and SIGLEC15+ macrophages. SIGLEC15+ macrophages, found in the peritubular space surrounding the testicular cords, interact with mesenchymal and endothelial cells via extracellular matrix proteins and are transcriptomically similar to osteoclasts, macrophages involved in bone reabsorption. Our results echo previous research showing that osteoclasts contribute to bone angiogenesis through a mechanism involving metalloproteinase-9 86 , an 580 extracellular matrix protein that is specifically expressed by SIGLEC15+ testicular macrophages. The spatial location and predicted interactions of these SIGLEC15+ macrophages suggest a possible role in mesonephric endothelial cell migration 77 . Indeed, SIGLEC15+ macrophages peak around 11-12 PCW, coinciding with the period of testicular endothelial remodelling. 585 The other macrophage population, microglia-like cells, are located primarily inside the testicular cords in close contact with Sertoli and germ cells and express immunomodulatory molecules. The transcriptome of these microglia-like cells is similar to those of fetal microglia and a subset of macrophages in the fetal skin. This transcriptional programme shared across 590 diverse tissues suggests that they have a protective function in dampening damaging immune responses. The testes need to maintain an immunomodulatory environment to prevent immune cells from attacking neoantigens appearing with the onset of spermatogenesis, which only starts during puberty 87,88 . Similarly, the anagen hair follicle in adults is considered a relatively immunosuppressive environment where antigen-presentation responses are strongly 595 downregulated 89 . This is the first description of immune cells derived from the embryo inside the human fetal cords. Future work will determine if these macrophages play a role protecting germ cells in postpubertal samples.
In females, the developing ovaries undergo dynamic cellular and morphological changes to 600 form a highly organised tissue structure that contains the primordial follicle reserve required for sustained female fertility. To systematically quantify the dynamic cell-cell interactions in the distinct spatial microenvironments, and to improve the accuracy of our predictions, we have developed a novel resource named CellSign, which links receptors to their downstream transcription factors. This module is now incorporated into CellPhoneDB.org v4.0. 605 After their arrival at the gonadal ridge, PGCs remain in the outer ovarian cortex, which is rich in extracellular matrix proteins and chemoattractants. Their pattern of adhesion molecules then alters, allowing them to migrate towards the outer medulla, where they will further differentiate. We show that pre-granulosa cells form a gradient in the ovarian tissue that will 610 regulate the distinct stages of germ cell development. As previously described in mice 39,44,45 , we define two waves of pre-granulosa cells. The first wave (preGC-I) derives from bipotent sLGR5 and will likely form the primary sex cords in females 90 . preGC-I cells produce estrogen, which is recognized by oogonia, indicated by the expression of estrogen receptors in oogonia STRA8+ (pre-meiotic). The second wave of LHX2+ pre-granulosa cells, derived from 615 the ovarian cortex, is distributed in the cortex (preGC-IIa-b) and transition zone, close to the medulla (preGC-IIc and developing granulosa). These will develop into the secondary sex cords. Cortical preGC-II (preGC-IIa-b) express CYP26B1, which inactivates retinoic acid, and ligands whose cognate receptors are expressed in PGCs. In contrast, transition zone preGC-IIc express ligands involved in initiation of meiosis (eg. retinoic acid by ALDH1A2) and oogenesis 620 (eg. BMP2, PTHLH). Markers of folliculogenesis (eg. NOTCH2/3) are only active in developing granulosa cells at 17 PCW or older. Human pre-granulosa subsets show a different transcriptomic profile compared to their murine counterparts, reinforcing the importance of continued study of their regulatory networks. Interestingly, we see functional conservation in estrogen and BMP2 production between murine and human pre-granulosa subsets, indicating 625 that the spatiotemporal functional compartmentalisation in the ovary is preserved across primates and rodents, even if the cell states and molecular regulators are different.
In conclusion, we have defined a comprehensive map of cells, molecular and cellular programmes of gonadal differentiation as a unique resource to understand gonadal function and 630 pathology. DSD are rare conditions associated with infertility and cancer. Many of these disorders involve pathology in other organs, (eg. kidney, heart, peripheral and central nervous systems) 91 . Our resource will help with prompt diagnosis of DSD by providing an accurate interpretation of sequencing datasets, finding novel therapeutic candidates, and improving clinical outcomes and quality of life. Our study will also help understand the pathogenesis of 635 reproductive cancers; for example, ovarian cancer recapitulates fetal programmes and germcell tumours originate during fetal development 92 . Finally, our atlas can be used as a blueprint for the design of protocols for in vitro gametogenesis from reprogrammed somatic cells, with relevance for infertility treatments. 640  proportional to scaled log-transformed expression. "o" = TFs whose binding motifs are 720 differentially accessible (i.e. TF can bind their potential targets), "a" = TFs whose targets are also differentially expressed (i.e. differentially activated TFs), and asterisk (*) = TFs that meet both "o" and "a'' conditions. f, Enrichment z-scores for peaks contained within CCANs including PCOS-associated single nucleotide polymorphisms (SNPs) with respect to the cell states identified in scATAC-seq data of human supporting cells. CCANs that contain multiple 725

Figure legends
SNPs are marked with a number of asterisks (*) that corresponds to the number of SNPs. g, Co-accessibility and coverage plot of the genomic region of the CCAN including the SNP rs8043701-A and gene TOX3. Viewpoint in the co-accessibility plot is set to the genomic coordinates of the SNP rs8043701-A. h, Dot plot showing the variance-scaled, log-transformed expression of TOX3 characteristic of pre-granulosa cell states. coelEpi = Coelomic epithelium 730 sLGR5 = supporting LGR5; sPAX8 = supporting PAX8; sKITLG = supporting KITLG; preGC = pre-granulosa cells; PV = Perivascular cells; DEG = differentially expressed gene.

Assignment of fetal developmental stage
Embryos up to eight post-conception weeks were staged using the Carnegie staging method 93 . At fetal stages beyond eight post-conception weeks, age was estimated from measurements of foot length and heel-to-knee length and compared with the standard growth chart 94 . A piece of skin, or where this was not possible, chorionic villi tissue was collected from every sample for 1000 Quantitative Fluorescence-Polymerase Chain Reaction analysis using markers for the sex chromosomes and the following autosomes 13,15,16,18,21,22, which are the most commonly seen chromosomal abnormalities. All samples were karyotypically normal.

Tissue processing
All tissues for sequencing and spatial work were collected in saline (HypoThermosol 1005 biopreservation media) and stored 4°C until processing. Tissue dissociation was conducted within 24 hours of tissue retrieval with the exception of tissues that were cryopreserved and stored in -80°C (see Supplementary Table 1).
We used the previous protocol optimised for gonadal dissociation 5 and available at protocols.io 95 . In short, tissues were cut into <1 mm 3 segments before being transferred to a 15ml tube. 1010 Tissues were digested with Trypsin/EDTA 0.25% between 5-15 minutes at 37°C with intermittent shaking. Digested tissue was passed through a 100 µm filter, and cells collected by centrifugation (500 x g for 5 minutes at 4°C). Cells were washed with flow buffer (PBS containing 5 % (v/v) FBS and 2 mM EDTA) prior to cell counting.

Single nuclei suspension
Single-nuclei suspension was produced from dissociated cells for doing scATAC-seq, following the manufacturers' instructions, and from frozen tissue sections for doing the multiomic snRNA-seq/scATAC-seq. For the latter, thick (300 µm) sections were 1025 cryosectioned, and kept in a tube on dry ice until subsequent processing. Nuclei were released via Dounce homogenisation as described in detail in the protocols.io 96 .

Tissue freezing
Fresh tissue samples of human fetal gonads were embedded in cold OCT medium and flash frozen using a dry ice-isopentane slurry. Protocol available at protocols.io 97 . 1035 Ovarian tissue collection from mouse embryos Developing ovaries were collected from E13.5 mouse embryos carrying the Oct4ΔPE-GFP transgene. Tissues were fixed in 4% (w/v) formaldehyde solution for 2 hours and at 4 °C. Samples were washed with PBS and afterwards sequentially incubated with 10% and 20% (w/v) sucrose at 4°C. After, samples were embedded in an optimal cutting temperature 1040 compound (OCT) and subsequently flash frozen using a dry ice-isopentane slurry. All experimental procedures were carried out in agreement with the project licence PE596D1FE issued by the UK Home Office and carried out in a Home Office designated facility.

Haematoxylin and Eosin (H&E) staining and imaging
Slides holding fresh frozen sections were removed from -80°C storage and let to air dry before 1045 being immersed in 10% Neutral Buffered Formalin for 5 minutes. After rinsing with RO water, slides were dipped in Mayer's Haematoxylin solution for 1.30 minutes. A series of submersions in RO water -with 4-5 changes of water and until water runs clear-was used to both rinse slides as well as blueing haematoxylin stain. 1% Aqueous Eosin was manually applied on top of the sections with a pipette and rinsed with RO water after 1-3 seconds. Slides were next 1050 immersed in a series of ethanol solutions (2 x 70% EtOH and 2 x 100% EtOH) to dehydrate tissue sections, before proceeding to 2 x Xylene baths. Slides were finally coverslipped and let to air dry before being imaged on Hamamatsu Nanozoomer 2.0HT digital slide scanner.

Multiplexed smFISH and high-resolution imaging
Large tissue section staining and fluorescent imaging was conducted largely as described 1055 previously 98 . Sections were cut from fresh frozen or fixed frozen samples embedded in OCT at a thickness of 10 µm using a cryostat, placed onto SuperFrost Plus slides (VWR) and stored at -80°C until stained. For FFPE samples, sections were cut at a thickness of 5 µm using a microtome, placed onto SuperFrost Plus slides (VWR), and left at 37°C overnight to dry and ensure adhesion. Tissue sections were then processed using a Leica BOND RX to automate 1060 staining with the RNAscope Multiplex Fluorescent Reagent Kit v2 Assay (Advanced Cell Diagnostics, Bio-Techne), according to the manufacturers' instructions. Probes may be found in Table S15. Prior to staining, human fresh frozen sections were post-fixed in 4% paraformaldehyde in PBS for 15 minutes at 4°C, then dehydrated through a series of 50%, 70%, 100%, and 100% ethanol, for 5 minutes each. Following manual pre-treatment, automated 1065 processing included epitope retrieval by protease digestion with Protease IV for 30 minutes prior to probe hybridisation. Mouse fixed frozen sections were subjected to the same manual pre-treatment described above. Subsequently, the automated processing for these sections included heat-induced epitope retrieval at 95°C for 5 minutes in buffer ER2 and digestion with Protease III for 15 minutes prior to probe hybridisation. 10x Genomics Visium library preparation and sequencing Ten micron cryosections were cut and placed in duplicate on Visium slides (beta product 1100 version). These were processed according to the manufacturer's instructions. Briefly, sections were fixed with cold methanol, stained with haematoxylin and eosin and imaged on a Hamamatsu NanoZoomer S60 before permeabilisation, reverse transcription and cDNA synthesis using a template-switching protocol. Second-strand cDNA was liberated from the slide and single-indexed libraries prepared using a 10x Genomics PCR-based protocol. 1105 Libraries were sequenced (1 per lane on a HiSeq4000), aiming for 300M raw reads per sample with read lengths 28cy R1, 8cy i7 index, 0cy i5 index, 91cy read 2.
Alignment and quantification of scRNA-seq data For each sequenced library, we performed read alignment to the 10x Genomics' GRCh38 3.1.0 reference genome, quantification and initial quality control (QC) using the Cell Ranger 1110 Software (version 3.1; 10X Genomics) using default parameters. Cell Ranger filtered count matrices were used for downstream analysis.
Downstream scRNA-seq analysis Quality filters, doublet detection, alignment of data across different batches, and clustering We used Scrublet for cell doublet calling on a per-library basis. We used a two-step diffusion 1115 doublet identification followed by Bonferroni-FDR correction and a significance threshold of 0.01, as described in 62 . Predicted doublets were not excluded from the initial analysis, but used afterwards to flag clusters with high doublet score.
We integrated the filtered count matrices from Cell Ranger and analysed them with Scanpy 1120 v.1.7.0, with the pipeline following their recommended standard practices. Briefly, we excluded genes expressed by less than three cells and excluded cells with fewer than 300 genes or with more than 20% mitochondrial reads. After converting the expression space to log(CPM/100 + 1), the object was transposed to gene space to identify cell cycling genes in a data-driven manner, as described in 62,74 . After performing PCA, neighbour identification and 1125 Louvain clustering, the members of the gene cluster including known cycling genes (CDK1, MKI67, CCNB2 and PCNA) were flagged as the data-derived cell cycling genes. Next, we identified the identifying highly variable genes, which were used for Principal Components Analysis (PCA). We corrected for the donor's effect using Harmony 32 on PCA space with theta equal to 0, keeping the remaining parameters default. Finally, we used the Harmony-corrected 1130 PCA matrix for neighbour identification, Louvain clustering and Uniform Manifold Approximation and Projection (UMAP) visualisation.

Annotation of scRNA-seq datasets
Identification, labeling and naming of the major cell types was carried by manual inspection of marker genes and interpretation of these based on the literature. We used the TF-IDF approach 1135 from the SoupX R package v.1.5.0 99 to identify genes specifically expressed in a cluster. Clusters with overall high doublet score or low counts number and no distinctive differentially expressed gene were flagged and discarded in further analysis.
Lineage-specific analysis 1140 Germ cells, supporting and mesothelial, mesenchymal, endothelial and epithelial cells were reanalysed independently using Seurat 3.2.2. For each lineage, we performed two analysis rounds. In both rounds, cells from a specific lineage (defined in the main analysis) were subsetted. Genes expressed by less than 3 cells were discarded. Counts were normalized, PCA was carried out on highly variable genes using Seurat defaults, and the donor effect was 1145 corrected by the Harmony tool (theta = 0). Corrected PCA dimension was used for neighbourhood identification, Louvain clustering and UMAP visualisation. We used TF-IDF 99 to identify cluster-specific genes. The first round of per-lineage analysis aimed to detect doublets and low QC cells not detected in the main analysis, and was characterized by high resolution clustering and identification of low quality cells. Briefly, we further annotated: 1) 1150 doublets, by identifying clusters with high scrublet score and expressing genes specific of multiple lineages; 2) low QC clusters, by identifying clusters with low counts number and without any distinctive differentially expressed genes; 3) low QC clusters with low counts number and high expression of genes specific of the erythroid lineage; and 4) cycling clusters composed mostly of cells in G2/M and S phase, as estimated with CellCycleScoring function 1155 in Seurat. The remaining cells (aka the clean per-lineage dataset) were re-analysed as described previously and annotated using known marker genes (when available) or the most distinctively expressed gene according to the TF-IDF approach.

Immune cells analysis 1160
Cell Ranger filtered count matrices of CD45+ enriched samples were processed using the Scanpy v.1.7.0 workflow described above for the main scRNA-seq analysis (Doublet detection, alignment of data across different batches, and clustering). These cells were then merged with the cluster of immune cells from the non-enriched samples. We corrected for donor effect using Harmony (theta = 0). We used the corrected PCA space to compute the neighborhood graph, 1165 Leiden clustering and UMAP visualisation. The resulting clustered manifold was preliminary annotated by transferring labels from a publicly available dataset of human fetal liver hematopoiesis 62 . Fetal liver raw sequencing data were downloaded from ArrayExpress (E-MTAB-7407), processed with Scanpy v.1.7.0 workflow described above for the main scRNAseq analysis and filtered based on the expression of CD45 (PTPRC) to exclude non-immune 1170 cells. We then trained a Support Vector Machine (SVM) model (sklearn.svm.SVC) on the filtered liver dataset and used it to predict cell types on our gonadal immune dataset. Predicted cell type annotations were validated or disproved by looking at the expression of known marker genes. 1175 To further characterize the novel gonadal macrophage populations, we performed functional enrichment analysis of Gene Ontology (Biological Process) terms with the Python interface of the g:Profiler 100 toolkit (gprofiler-official v.1.0.0) using the top 100 marker genes identified with TF-IDF as query and all the genes as background. Projection of fetal microglia cells 63 onto our immune dataset was done using a SVM model as well as scmap 101 following the 1180 tutorial (http://bioconductor.org/packages/release/bioc/vignettes/scmap/inst/doc/scmap.html).

Cross-species comparison
We compared our single-cell signatures of germ and supporting cells to those of other 1195 mammalian species by re-analysing two published datasets of mouse and one of macaque fetal gonads. The macaque dataset included fetal ovaries at stages E84 and E116 38 . One mouse dataset included fetal ovaries from E11.5 to P5 39 , while the other contained germ cells from both fetal ovaries and testes from E10 to E16 40 . The distinct datasets were downloaded from the Gene Expression Omnibus database (GSE136220, GSE136441, GSE149629). We re-1200 analysed the datasets using the Scanpy v.1.7.0 workflow, including batch correction with Harmony, as described above for the main scRNA-seq analysis. Macaque and mouse genes were converted to human genes using ENSEMBL Biomart multi species comparison filter. Main cell populations were annotated based on the expression of marker genes described by the authors. For the comparison of germ cells, we subsetted the datasets to keep germ cells 1205 only, and plotted our fetal human signatures onto them. For the comparison of ovarian supporting cells, we trained a SVM model (sklearn.svm.SVC) on human ovarian supporting cells and projected the data onto the macaque and mouse datasets. Additionally, we plotted human genes in the mouse and macaque dataset to validate our models, and found additional functional relationships between the development of gonadal supporting cells in mammals. 1210

Alignment, quantification, and quality control of ATAC data
We processed scATAC-seq libraries (read filtering, alignment, barcode counting, and cell calling) with 10x Genomics Cell Ranger ATAC pipeline (version 1.2.0) using the pre-built 10x's GRCh38 genome (version 3.1.0) as reference. We called the peaks using an in house 1215 implementation of the approach described in Cusanovich et al. 103 (available at https://github.com/cellgeni/cellatac, revision 21-099). In short, the genome was broken into 5 kb windows and then each cell barcode was scored for insertions in each window, generating a binary matrix of windows by cells. Matrices from all samples were concatenated into a unified matrix, which was filtered to retain only the top 200K most commonly used windows per 1220 sample. Using Signac (https://satijalab.org/signac/ version 0.2.5), the binary matrix was normalised with TF-IDF followed by a dimensionality reduction step using Singular Value Decomposition (SVD). Latent Semantic Indexing (LSI) were clipped at +/-1.5. The first LSI component was ignored as it usually correlates with sequencing depth (technical variation) rather than a biological variation 103 . The 2-30 top remaining components were used to perform 1225 graph-based Louvain clustering. Next, peaks were called separately on each cluster using macs2 104 . Finally, peaks from all clusters were merged into a master peak set (i.e. peaks overlapping in at least one base pair were aggregated) and used to generate a binary peak by cell-matrix, indicating any reads occurring in each peak for each cell. 1230

Downstream scATAC-seq analysis
Quality filters, alignment of data across different batches, and clustering To obtain a set of high quality peaks for downstream analysis, we filtered out peaks that (i) are included in the ENCODE blacklist, (ii) have a width outside the 210-1500bp range and (iii) are accessible in less than 4% of cells from a cellatac cluster. Low quality cells were also removed 1235 by setting to 5.5 the minimum threshold for log1p transformed total counts per cell.
We adopted the cisTopic approach 105,106 for the core of our downstream analysis. cisTopic employs Latent Dirichlet Allocation (LDA) to estimate the probability of a region belonging to a regulatory topic (region-topic distribution) and the contribution of a topic within each cell 1240 (topic-cell distribution). The topic-cell matrix was used for constructing the neighborhood graph, computing UMAP projections and clustering with the Leiden algorithm. Donor effects were corrected using Harmony (theta = 0), with the exception of the germ cells dataset which already exhibited good batch mixing in the uncorrected manifold. Cell doublets were identified and removed using scrublet 107 . 1245 Gene activity scores Next, we generated a denoised accessibility matrix (predictive distribution) by multiplying the topic-cell and region-topic distribution and used it to calculate gene activity scores. To integrate them with scRNA-seq data, gene activity scores were rounded and multiplied by a factor of 1250 10^7, as previously described 106 .
Cell type annotation To annotate cell types in scATAC-seq data, we first performed label transfer from scRNA-seq data of matched individuals. We used Canonical Correlation Analysis as dimensionality 1255 reduction method, vst as selection method, 3000 variable features and 25 dimensions for finding anchors between the two datasets and transferring the annotations 34 . The predicted cell type annotations by label transfer were validated by importing annotations of the combined snRNA-seq/snATAC-seq profiling data. 1260 Cell type-specific cis-regulatory networks Co-accessible peaks in the genome and cis-co-accessibility networks (CCANs) were estimated using the R package Cicero 36 v. 1.3.4.11 with default parameters. We then filtered the denoised accessibility matrix from cisTopic to keep only the peaks included in CCANs. The resulting matrix was further processed to average cells by cell type and peaks by CCAN. Finally, we z-1265 scored the matrix across CCANs and visualized the separation of CCANs by cell type by hierarchical clustering and plotting the heatmap.
GWAS enrichment with single-cell annotations We downloaded the list of polycystic ovary syndrome associated SNPs from the GWAS 1270 Catalog (Experimental Factor Ontology id = EFO_0000660) and filtered it to keep only the variants exceeding genome-wide significance (p-value = 5 x 10^8). Bed files were created from both the SNPs and the CCANs, and we used bedtools 108 to intersect the two sets of bed files. Z-scores of CCANs including a PCOS-associated SNP were visualized with a heatmap to assess which cell types are most likely to be affected by these SNPs. Moreover, we specifically 1275 focused on the CCAN containing the variant rs8043701-A and investigated how this variant could disrupt the co-accessibility pattern of the CCAN (using the plot_connections function in Cicero and setting as viewpoint the genomic coordinates of variant rs8043701-A).

Alignment, quantification, and quality control of Visium data
For each 10X Genomics Visium sequencing data, we used Space Ranger Software Suite 1280 (v.1.2.1) to align to the GRCh38 human reference genome (official Cell Ranger reference, version 2020-A) and quantify gene counts. Spots were automatically aligned to the paired H&E images by Space Ranger software. All spots under tissue detected by Space Ranger were included in downstream analysis.

Downstream analysis of 10x Genomics Visium data
1285

Location of cell types in Visium data
To spatially locate the cell states on the Visium transcriptomics slides, we used the cell2location tool v.0.3 33 . As reference, we used scRNAseq data from individuals of the same sex and gestational stage. We used general cell annotations from the main analysis, with the exception of the major gonadal lineages (germ, supporting and mesenchymal) for which we 1290 considered the identified sub-populations. We used default parameters with the exception of cells_per_spot which was set to 20. Each Visium section was analysed separately. Results were visualized following the cell2location tutorial.

1295
We updated the CellphoneDB database to include: 1) more manually curated protein cell-cell interactions (n = 1,852 interactions) and 2) cell-cell interactions involving non-protein ligands such as steroid hormones and other small molecules (n = 194). For the latter, we used the last bona fide enzyme in the synthesis pathway. Additionally, we have added a new module called CellSign which links receptors in CellphoneDB to their known downstream transcription 1300 factor, which were manually curated and the relevant pubmed reference number recorded. The code is available at https://github.com/ventolab/Cellphonedbv3, and a Python package will be provided in the revised version.
To retrieve interactions between supporting and other cell populations identified in our gonadal 1305 samples, we used an updated version of our CellPhoneDB 76,109 approach described in 27 . In short, we retrieved the interacting pairs of ligands and receptors meeting the following requirements: 1) all the protein members were expressed in at least 10% of the cell type under consideration; and 2) at least one of the protein members in the ligand or the receptor is a differentially expressed gene, with an adjusted p-value below 0.01 and a log2 fold change 1310 above 0.2. To account for the distinct spatial location of cells, we further classified the cells according to their location in the fetal ovaries (cortex, medulla) as observed by Visium and smFISH. We filtered cell-cell interactions to exclude cell pairs that do not share the same location. 1315

Transcription Factor (TF) analysis
To prioritize the TFs relevant for a cell state in a lineage, we integrate three measurements: (i) expression levels of the TF and (ii) the activity status of the TF measured from (iia) the expression levels of their targets (described below in "Transcription factor activities derived from scRNAseq") and/or (iib) the chromatin accessibility of their binding motifs (described 1320 below in "Transcription factor motif activity analysis from scATACseq"). Plots include TFs those meeting the following criteria: 1) TF is differentially expressed (i), with log2 fold change greater than 0.5 and adjusted p-value < 0.05 and 2) TF is differentially active (ii), with log2 fold change greater than 1 and adjusted p-value < 0.05 in at least one of the TF activity measurements (iia/iib). 1325 Transcription factor differential expression (from scRNAseq) We compute differential expression using the Logistic Regression approach implemented in the FindAllMarkers function in Seurat, in a one-vs-all fashion. 1330 Transcription factor activities derived from scRNAseq We estimated protein-level activity for human Transcription factors (TF) as a proxy of the combined expression levels of their targets. Target genes were retrieved from Dorothea 110 , an orthogonal collection of TF targets compiled from a range of different sources. Next, we estimated TF activities for each cell using Viper 111 , a GSEA-like approach, as implemented in 1335 the Dorothea R package and tutorial 112 . Finally, to identify TFs whose activity is upregulated in a specific cell type, we applied Logistic Regression from Seurat onto the z-transformed cell x TF activity matrix in a one-vs-all fashion. TFs with adjusted p-value below 0.05 and average log2 fold change above 1 in a cell type were considered as differentially active . 1340 Transcription factor motif activity analysis from scATACseq Transcription factor motif activities were computed using chromVar 35 v. 1.12.2 with positional weight matrices from JASPAR2018 113 , HOCOMOCOv10 114 , SwissRegulon 115 , HOMER 116 . chromVar returns a matrix with binding activity estimates of each TF in each cell, which we used to test for differential TF binding activity between cell types in a one-vs-all fashion with 1345 Logistic Regression (FindAllMarkers function in Seurat).