Introduction

Early mammalian embryonic gonads consist of primordial germ cells (PGCs) and surrounding gonadal somatic cells. PGCs are specified in the epiblast near the proximal allantois. In humans, they begin migrating at approximately 4 weeks of development¹ (embryonic day E7.5 in mice², E20 in bovine³) and travel through the hindgut endoderm before colonizing the genital ridge. At this stage, the genital ridge is bipotential and has the capacity to develop into either testes or ovaries⁴. Gonadal somatic cells arise from the genital ridge at around 32 days of pregnancy in humans (E10.5 in mice, E27-E31 in bovine), in which coelomic epithelial cells differentiate into two major somatic precursor lineages, including supporting cell precursors and steroidogenic cell precursors⁵.

In males, upregulation of the SRY-SOX9-FGF9 axis initiates male sex determination and drives somatic cells to differentiate into pre-Sertoli cells⁶. In females, in the absence of SRY, the RSPO1/WNT4/β-catenin and FOXL2 pathways are activated to promote gonadal somatic cell differentiation into pre-granulosa cells⁷. As pre-granulosa cells differentiate, PGCs initiate meiosis and eventually arrest at the diplotene stage of prophase I before puberty⁸. In parallel, PGCs undergo extensive epigenetic reprogramming, characterized by global DNA demethylation and chromatin remodeling that establishes a permissive transcriptional landscape⁹.

PGC proliferation and early germ cell differentiation in females are controlled by conserved regulators such as DAZL^10,11 and DDX4^12,13, which are essential for germ cell maintenance and the acquisition of mitotic competence across mammals. While core features of gonadal development are broadly conserved across mammals, important species-specific regulatory mechanisms have also been reported. Meiosis is subsequently initiated by STRA8, a factor essential for meiotic entry in mice¹⁴, whereas in humans, meiotic initiation is associated with STRA8 expression but appears to involve a more complex regulatory network, with partial dependence on retinoic acid¹⁵. Despite these differences, core meiotic regulators such as SYCP3^13,14 and REC8^16,17 are conserved across mammals.

Early gonadal development has important implications for both human health and livestock production. Impairments in these processes can result in disorders of sex development (DSDs) or infertility in mammals¹⁸. For example, mutations in SRY can cause 46,XY gonadal dysgenesis and sex reversal in humans¹⁹. In livestock species, abnormalities in gonadal development, such as ovarian dysgenesis or ovarian hypoplasia²⁰, represent the most common cause of reproductive failure. In livestock production, early gonadal development has long attracted research interest due to its direct impact on puberty onset and its potential to improve reproductive efficiency. However, limited mechanistic understanding has restricted its practical application.

Recent advances in single-cell RNA sequencing (scRNA-seq) have enabled high-resolution characterization of gonadal development in humans²¹, mice²², and livestock²³, revealing both conserved and species-specific transcriptional dynamics in germ cells and somatic lineages during sex determination and ovarian differentiation. By resolving cellular heterogeneity, scRNA-seq facilitates the identification of rare or transient cell populations at a limited developmental window, such as supporting cell intermediates²¹, and the reconstruction of lineage-specific developmental trajectories for both germ cells and somatic cells. Despite these advances, cross-species comparison of gonadal development remains technically and experimentally challenging. For example, each species has differences in developmental timing, which makes it hard to align collected samples for equivalent stages. A limited number of one-to-one orthologous genes map across multiple species, which further challenges the number of genes that can be compared. Moreover, technical variability in single-cell data generation and processing from different studies further hinders integrative analyses. As a result, despite scRNA-seq having been generated in developing gonads of multiple mammalian species, a robust and systematic framework for delineating conserved and divergent regulatory mechanisms of gonadal development across mammals is still lacking.

In this study, we generated a bovine RNA-seq dataset covering six time points during early ovarian development (E38–E112) and performed a cross-species comparative scRNA-seq analysis with stage-matched human (PCW6–16) and mouse (E11.5–E18.5) datasets. Our objectives were: 1) compare the dynamics of germ and somatic cell populations across species during major events in early development, 2) identify conserved and species-specific regulators involved in germ cells and granulosa differentiation, and 3) characterize cell-cell communication between germ cells and granulosa cells during fetal ovarian development. Together, this analysis offers a comparative framework to better understand how the early female gonad is regulated at the molecular level during development across mammals, revealing conserved and species-specific regulatory programs and cellular dynamics.

Methods

Ethics Statement

All animal procedures were conducted in accordance with the National Research Council’s Guide for the Care and Use of Laboratory Animals and were approved by the Institutional Animal Care and Use Committee (IACUC) of Northwest A&F University (protocol IACUC2024-1204, approved on 10 December 2024). Fetal tissues were obtained from a local slaughterhouse as previously described²⁴.

Fetal Bovine Gonad Collection

Estrus in donor cows was monitored daily to determine the timing of pregnancy. Artificial insemination was performed once at the onset of estrus. Uteri were collected from cows maintained under natural conditions and transported to a sterile laboratory for processing. Embryos were dissected, and fetal gonads were isolated under sterile conditions. Morphologically intact fetal gonads with sufficient cell viability were collected from female embryos at developmental stages E38, E46, E73, E83, E92, and E112.

Single-Cell Suspension Preparation of Fetal Gonads

To isolate germ cells and gonadal somatic cells, fetal gonads were dissected under a stereomicroscope and placed into 1.5-mL tubes containing an enzymatic digestion solution (1 mg/mL Collagenase IV, 1 mg/mL Hyaluronidase, 0.25% Trypsin, and 1 mg/mL DNase I). Samples were incubated at 37°C for 30 minutes, after which serum was added to stop the enzymatic reaction. The resulting cell mixture was gently pipetted, passed through a 40-μm strainer, and resuspended in 0.1% serum/PBS to generate single-cell suspensions for 10x Genomics scRNA-seq.

Single-Cell RNA Sequencing

Single-cell suspensions with >80% viability were adjusted to a final concentration of 700–1200 cells/μL. For each sample, single-cell encapsulation and library preparation were performed at Novogene Bioinformatics Technology Co., Ltd. (Tianjin, China) using the MobiCube® High Throughput Single Cell 3’ RNA-Seq Kit (v2.1), following the manufacturer’s instructions. Cells or nuclei were loaded into the microfluidic chip of the Chip A Single Cell Kit (MobiDrop, cat. no. S050100201) and processed on the MobiNova-100 system (MobiDrop, cat. no. A1A40001) to generate droplets. Each droplet contained a single cell or nucleus together with a gel bead carrying millions of oligonucleotides with unique cell barcodes. After droplet formation, light-mediated cleavage using the MobiNovaSP-100 (MobiDrop, cat. no. A2A40001) released the oligos into the reaction mix. mRNAs were captured by oligo(dT) within the gel beads, followed by reverse transcription and cDNA amplification. Libraries were constructed using the High Throughput Single Cell 3’ RNA-Seq Kit v2.0 (MobiDrop (Zhejiang) Co., Ltd., cat. no. S050200201) and the 3’ Single Index Kit (MobiDrop (Zhejiang) Co., Ltd., cat. no. S050300201). Sequencing was performed on an Illumina NovaSeq platform using 150-bp paired-end reads at Novogene (Tianjin, China).

Single-Cell RNA Sequencing Data Processing

Raw sequencing data were processed using Cell Ranger (v3.1, 10x Genomics)²⁵. For bovine samples generated in this study, sequencing reads were aligned to the bovine reference genome (ARS-UCD2.0) using Cell Ranger according to the manufacturer’s instructions. Human (GRCh38) and mouse (GRCm38) reference genomes were used for processing the corresponding datasets. Species-specific reference genomes were built using the Cell Ranger mkref pipeline. The filtered feature–barcode matrices generated by Cell Ranger were used for downstream analyses.

Quality Control and Cell Type Annotation

The quality control was performed using the Seurat package (v5.2.1)²⁶. Cells expressing fewer than 200 genes and genes detected in fewer than 3 cells were excluded. Cells with nFeature_RNA or nCount_RNA values outside ±3 median absolute deviations (MADs) from the median were classified as outliers and removed. Additionally, cells with mitochondrial gene expression exceeding 10% were excluded. Doublets were identified and removed using scDblFinder (v1.61.0)²⁷. One-to-one orthologous genes from Ensembl were used to integrate single-cell RNA-seq datasets from bovine, human, and mouse. Data were normalized using Seurat, and the top 5,000 highly variable genes were selected for downstream analyses. During data normalization, sequencing batches, cell cycle scores, nFeature_RNA, and the percentage of mitochondrial and ribosomal genes were regressed out. For identifying cell clusters, the optimal resolution was selected with the aid of the clustree (v0.5.1) package²⁸ and the dimensionality reduction was performed using the top 25 principal components. To annotate cell clusters, the FindAllMarkers function was used with default parameters to identify cluster-specific marker genes. Each cluster was manually annotated using canonical markers to identify mammalian ovarian cell types.

Pseudotime Trajectory Analysis

To reconstruct developmental trajectories of germ cells and granulosa cells, pseudotime analysis was performed using Monocle3 (v1.3.7)²⁹ following the standard workflow. Analyses were conducted separately for each species. For each cell type, cells enriched for early-stage markers were selected to define the root of the trajectory, and pseudotime values were inferred using Monocle3. Trajectory-associated genes were identified by spatial differential expression analysis using the graph_test() function (q value < 0.05). For each cell type, trajectory-associated genes were grouped into 3-5 modules along pseudotime by unsupervised hierarchical clustering and visualized using ClusterGVis (v0.1.4)³⁰. Gene sets identified from each species were then compared using Venn diagrams to determine shared and species-specific trajectory-associated genes.

Gene Ontology (GO) enrichment analysis

Trajectory-associated genes identified from the pseudotime analysis were subjected to Gene Ontology (GO) enrichment analysis using ShinyGO (v0.85.1)³¹. GO analysis was performed separately for each species.

Transcriptional Regulatory Network Analysis

Gene regulatory network analysis was performed using pySCENIC (v0.12.1)³² for germ cells and granulosa cells, with analyses conducted separately for each species. Raw count matrices were used to infer transcription factor–centered gene co-expression networks using GRNBoost2. These co-expression modules were subsequently filtered by cis-regulatory motif enrichment analysis to define regulons comprising transcription factors and their predicted target genes. Regulon activity was quantified at the single-cell level using AUCell, generating regulon activity (AUC) scores. Transcription factors exhibiting conserved or species-specific regulon activity were identified during germ cell and granulosa cell developmental progression.

Cell-Cell Communication Analysis

Cell-cell communication analysis was performed using the CellChat package (v2.1.0)³³. Cell–cell interactions mediated by ligand–receptor pairs were inferred with the default parameters. The normalized counts were used as input.

SVM-Based Cross-Species Cell Type Comparison

To compare cell-type transcriptomic similarity across species, a linear support vector machine (SVM)³⁴ was trained on the human dataset and applied to mouse and bovine cells. To reduce class imbalance, human cell types were downsampled to match the smallest cell-type population, 305 cells. 75% of the data were used for training the model and 25% of the data were used to test the model. This procedure was repeated using bootstrapping to ensure the robustness of model performance. Raw UMI counts were normalized to counts per million (CPM) using edgeR³⁵ and were log1p-transformed. The top 2,000 highly variable genes were selected based on variance in the human dataset and used as features for model training and projection. The SVM was trained using LiblineaR³⁶ with annotated human cell types as labels and subsequently used to project mouse and bovine cells. For each cell, the predicted probability corresponding to its annotated cell type was extracted and summarized by cell type to assess relative cross-species similarity.

To further examine species-specific transcriptional features within granulosa cells, an additional multiclass SVM was trained on granulosa cells from human, mouse, and bovine. Species-associated genes were identified by ranking the absolute values of the learned linear coefficients for each class. In addition, to infer the identity of the unclassified bovine-specific cell population, an SVM was trained on bovine cells excluding this population and used to predict the most similar bovine cell types.

Results

Integration of single-cell RNA-seq data from early female gonads across bovine, human, and mouse

To characterize the cellular landscape of early female gonadal development in cattle, we generated single-cell RNA sequencing (scRNA-seq) profiles from fetal bovine ovaries across six gestational stages (Figure 1A). These stages encompass key developmental milestones: sex determination (E38, E46), pre-granulosa cell differentiation (E73, E83, E92, and E112), and development of early germ cells, including primordial germ cells (PGCs), meiotic germ cells and early oocytes. Here, we collectively refer to germ cells across all developmental stages, including PGCs, meiotic germ cells, and early oocytes as “germ cells”. To enable cross-species comparison, we integrated publicly available scRNA-seq datasets from human fetal ovaries (post-conception weeks PCW6, PCW7, PCW10, PCW12, PCW13, and PCW16)²¹ and mouse ovaries (E11.5, E12.5, E14.5, E16.5, and E18.5)²², selected to represent developmentally comparable time windows matched to the bovine samples (Figure 1A).

Figure 1.

Following integration based on 16,073 one-to-one orthologous genes, the combined dataset comprised 140,441 single cells across all three species (Supplementary Table 1), including 52,928 bovine, 48,874 human, and 38,639 mouse cells. After quality control, 107,930 high-quality cells were retained for downstream analysis (Supplementary Table 1). Unsupervised clustering followed by manual annotation guided by marker genes established in human and mouse gonadal development²¹ identified 11 major cell types (Figure 1B, Supplementary Table 2), including germ cells, granulosa cells, mesenchymal cells, mesothelial cells, epithelial cells, erythroid cells, endothelial cells, immune cells, perivascular (PV) cells, glomerular mesonephros cells (transient embryonic kidney), neural cells, as well as a population of unclassified bovine-specific cells. For example, germ cells were identified by high expression of DAZL and DDX4 (Figure 1C), granulosa cells by KITLG and AMHR2 (Figure 1C), epithelial cells by PAX8, and mesenchymal cells by DCN (Supplementary Figure 1).

Across all species and developmental stages, mesenchymal cells, and granulosa cells together constituted the majority of the cellular population in the early development, highlighting their central roles during early ovarian development (Figure 1B and 1D). Interestingly, female germ cells became more predominant after the initiation of the meiosis stage, especially in bovine (E73) and mouse (E14.5) datasets. Moreover, the unclassified bovine-specific cell population was detected at all sampled stages, with a higher relative abundance at early gestational time points (E38, E46, and E73), suggesting the presence of species-specific cellular states during early bovine gonadal development (Figure 1B and 1D).

Dynamics of epigenetic regulators during early ovarian development

To investigate epigenetic reprogramming during early development of female germ cells and granulosa cells, we analyzed the gene expression dynamics of key regulators associated with DNA methylation (DNAme), histone modifications, and chromatin remodeling across the 3 species (Figure 2; Supplementary Figure 2).

Figure 2.

In female germ cells (Figure 2), for DNAme regulators (Figure 2A), we observed high expression of the maintenance DNA methyltransferase DNMT1³⁷ and upregulation of the DNA demethylation factor TET1³⁸ at early developmental stages (E46 in bovine; PCW7 in human and E12.5 in mouse) followed by decreased expression of DNMT1 and TET1 at later stages (Figure 2A). This pattern is consistent with the well-characterized wave of global DNA demethylation followed by remethylation during germ cell development³⁹.

Histone modification-related regulators also exhibited dynamic changes during early germ cell development (Figure 2B). EZH2, a core component of the Polycomb repressive complex 2 (PRC2) responsible for H3K27 trimethylation⁴⁰, showed sustained but gradually decreasing expression in mouse female germ cells (E11.5-E18.5). Histone demethylases KDM3A⁴¹ exhibited continuous activation in mouse germ cells after sex determination (E12.5-E18.5). Another histone demethylases KDM3B⁴¹ showed transient activation in bovine (E46) and mouse (E11.5). SMARCA1, a chromatin remodeler⁴², exhibited transient upregulation at early stages in bovine (E46) and human (PCW7), but not in mouse. In addition, SUZ12, another core component of the PRC2 complex required for its stability and catalytic activity⁴³, was upregulated in bovine germ cells (E46-E83) but showed decreasing expression in mouse female germ cells (E12.5-E18.5).

In granulosa cells, for DNAme regulators (Supplementary Figure 2A), expression of DNMT1³⁷ showed a progressive decrease following sex determination in bovine (E46 onward) and human (PCW7 onward). In contrast, mouse DNMT1 showed transient activation at E12.5 and E16.5. HELLS, a chromatin remodeler involved in DNA methylation maintenance and heterochromatin organization⁴⁴, showed transient activation at E46 in bovine and at E11.5, E12.5 and E16.5 in mouse. In contrast, TET1 displayed a gradual increase in expression after E73 in bovine, whereas in human and mouse it showed only transient upregulation at PCW13 and E16.5, respectively. UHRF1, another DNA methylation maintenance factor⁴⁵, showed transient activation at E11.5 in mouse. Together, these patterns suggest a shift from DNA methylation maintenance toward increased demethylation activity during granulosa cell development, although this shift varies across species in both timing and extent.

Histone modification-related regulators such as EZH2 showed dynamics similar to DNMT1 in each species, indicating the consistency in methylome erasure and reestablishment during development (Supplementary Figure 2B). In contrast, KDM3A and KDM6B showed distinct species-specific dynamics. In bovine, both genes were highly expressed at early stages (E38 and E46) and subsequently declined, whereas in human their expression was high at early development, showed a transient decrease after PCW12 and increased afterwards. In mouse, both KDM3A and KDM6B showed a sustained increase from E12.5 to E18.5. A similar divergence was observed for SMARCA1, which was highly expressed at early stages and gradually decreased in bovine, remained consistently high in human, and showed an increased expression pattern in mouse. Moreover, SUZ12 showed broader activation windows across species, which sustained high expression in bovine, a decrease then increase dynamic pattern in human, and a pronounced peaked expression at E16.5 in mouse. Together, these results suggest that epigenetic changes may also occur in granulosa cells during early development, although it has been more extensively characterized in germ cells.

Conserved and species-specific transcriptional dynamics during early germ cell development

To reconstruct the developmental progression of early germ cells, we performed pseudotime trajectory analysis of germ cells for each species (Figure 3A; Supplementary Figure 3). Across all three species, germ cell development followed a consistent transition from pluripotency-associated states toward meiotic initiation and differentiation. A set of genes has been identified that exhibited conserved dynamic expression patterns during this process (Figure 3A-B). For example, early development was marked by high expression of pluripotency and germline specification factors, including POU5F1⁴⁶, NANOG⁴⁷, TFAP2C⁴⁸, and SOX17⁴⁹, consistent with an undifferentiated germ cell identity. As the development progressed, these markers were gradually downregulated, coinciding with the activation of germ cell maturation regulators such as KIT⁵⁰, and RNA-binding factors ELAVL2⁵¹. At later stages, meiotic entry and germ cell differentiation were characterized by the upregulation of meiotic genes, including SYCP1⁵², SYCP2⁵³, SYCP3⁵⁴, TEX30⁵⁵, ZCWPW1⁵⁶, and SMC1B⁵⁷ (Figure 3A; Supplementary Figure 3B). We then overlapped trajectory-associated genes across the three species and identified a total of 34 conserved, highly dynamic genes, including NANOG, POU5F1, SYCP1, and JUN (Figure 3B). Beyond these conserved regulators, genes involved in cell-cycle control such as UBE2C⁵⁸ and TOP2A⁵⁹, translational regulation such as CPEB1⁶⁰, and broader differentiation and signaling processes represented by JUN⁶¹ and NR6A1⁶², also displayed dynamic expression changes along pseudotime across all three species. In addition, a substantial number of genes displayed species-specific expression dynamics, including 162 bovine-specific, 36 human-specific, and 566 mouse-specific genes (Figure 3B). These differences indicate species-specific divergence in gene regulatory features despite a conserved core trajectory of germ cell development. For example, GO analysis of trajectory-associated genes in each species revealed species-specific pathway enrichment, with bovine female germ cells enriched for pathways related to cell morphogenesis, reproductive process, and neurodevelopment, human female germ cells enriched for cell division, mitotic cell cycle, and DNA replication processes, and mouse female germ cells enriched for meiosis-related pathways, including homologous chromosome pairing and segregation and meiotic cell cycle (Figure 3C). These differences likely reflect variation in developmental timing and sampling across species.

Figure 3.

To further characterize regulatory programs during germ cell differentiation, we defined transcription factor regulon activity in each species (Figure 3D). This analysis identified a set of conserved active regulons across species, including NANOG, WT1, GATA4, GATA6, PBX1, MAF, and HMGA2. Several of these regulators have been previously implicated in male gonadal or germline-related contexts. For example, PBX1, which shows the strongest signal in cattle, has been reported as a master regulator of Leydig cell differentiation and steroidogenesis-related function in mice⁶³. Other regulators include WT1, which is required for testis development⁶⁴, MAF, which contributes to germline stem cell maintenance⁶⁵, and HMGA2, which has been associated with postpubertal testicular germ cell tumour⁶⁶. How these regulators contribute to early female gonadal development remains to be determined.

In addition, several regulons exhibited species-specific activity patterns. For example, we found the enriched ELF1, ELK3, and RFX2 activity in bovine germ cells, elevated MAFF and TEAD4 activity in mouse germ cells, and increased NFE2L3, STAT1, and E2F1 activity in human germ cells. RFX2 is also a key transcriptional regulator of mouse spermiogenesis⁶⁷, although its potential role in early female germ cell development in bovine remains unclear. STAT1 is likely to be responsible for initiating transcription network in human fetal germ cells⁶⁸. E2F1 has been shown to regulate germ cell apoptosis during the first wave of spermatogenesis in mouse⁶⁹. Together, these results demonstrate that while mammalian germ cell development follows a broadly conserved trajectory across species, it is accompanied by both shared and species-specific gene expression programs and transcription factor regulons that reflect divergence in gene regulatory control.

Conserved and species-specific transcriptional dynamics during early granulosa cell development

Similarly, we performed pseudotime trajectory analysis for granulosa cells to study their development separately for each species (Figure 4A; Supplementary Figure 4). In bovine granulosa cells, early development was marked by WNT signaling pathway modulator SFRP⁷⁰, cell-cell adhesion regulator CDH4⁷¹, and DCN, which has a structural role in ovarian extracellular matrix⁷². At later stages, granulosa cells showed activation of genes associated with endocrine maturation and follicular support, including INHA⁷³, CYP11A1⁷⁴, and PAPPA⁷⁵. In human granulosa cells, early to intermediate stages showed expression of GPC3, which encodes a glypican found in the extracellular matrix²², and MDK, a growth factor that is thought to support postnatal follicle growth in human ovaries⁷⁶. Differentiated granulosa cells were marked by activation of key granulosa cell markers, including FOXL2⁷⁷ and GJA1⁷⁸, alongside IGF pathway modulation through IGFBP2⁷⁹ and retinoic acid metabolism via RDH10⁸⁰. In mouse granulosa cells, early developmental stage was also characterized by expression of developmental and growth factor associated regulators such as GPC3, MDK, and LHX9, a LIM homeobox transcription factor implicated in gonadal development and ovarian function in mouse models⁸¹. Late developmental stage in mice showed expression of genes involved in key signaling pathways, such as DHH, a Hedgehog ligand expressed in granulosa cells with roles in follicle signaling and theca differentiation⁸², and granulosa cell maturation was marked by upregulation of FOXL2, similar to humans (Figure 4A).

Figure 4.

Some of these genes were conserved across developmental stages in all three species, including FOS and JUNB, both members of the AP-1 transcription factor family⁸³. These factors are induced in granulosa cells in response to gonadotropins and contribute to follicle maturation and ovulatory gene regulation in mice⁸⁴ (Figure 4B). Moreover, species-specific trajectory-associated genes were also observed with specific GO term enrichment, further indicating their regulatory divergence. Specifically, we observed that bovine granulosa cells were enriched for tube morphogenesis and neurogenesis, while human cells prioritized cytoplasmic translation, and mouse cells emphasized muscle tissue development (Figure 4C).

Transcription factor regulon analysis also identified FOS and JUNB as the most enriched common active regulons across all three species (Figure 4D), especially in humans. In contrast, several regulons exhibited species-specific activity patterns among granulosa cells. In bovine, HLTF, a DNA repair-associated factor⁸⁵, and PBX1, a developmental transcriptional co-factor⁸⁶, showed elevated activity. In mice, JUND⁸³, a member of the AP-1 family, together with FOS and JUNB, exhibits increased regulon activity.

Cell-cell interaction features between germ cells and gonadal somatic cells

To characterize how cell-cell communication among female gonadal somatic cell populations and early germ cells contributes to the gonadal differentiation and gametogenesis, we defined ligand-receptor interactions among cell types within each species (Figure 5). Overall, interactions involving germ cells were relatively weak compared with those among gonadal somatic cell populations, which exhibited more frequent intercellular communications, a pattern consistently observed across all three species (Figure 5A).

Figure 5.

To further dissect germ cell-associated signaling, we first designated germ cells as the signaling source and examined ligand-receptor interactions between them and other gonadal cell types (Figure 5B). Several ligand-receptor pairs were conserved across species. For instance, MDK-SDC4 interactions, involving the growth factor Midkine (MDK)⁸⁷ and its co-receptor Syndecan-4 (SDC4)⁸⁸, were detected between germ cells and epithelial cells in all three species and were additionally shared between germ cells and granulosa cells in both human and mouse. Moreover, NECTIN3-NECTIN2 interactions, which regulate cell-cell adhesion and junction formation⁸⁹, were conserved between cattle and mouse and occurred between germ cells and multiple somatic cell types, including endothelial cells, granulosa cells, mesothelial cells, and perivascular cells (PV).

Species-specific ligand-receptor interaction pairs were also identified, for example, COL1A2-SDC4 and COLA1- SDC4 were only observed to signal out from human germ cells to either epithelial cells, granulosa cells, or mesothelial cells. EFNA5-EPHA5 were only sent from bovine germ cells into granulosa cells or mesothelial cells. FN1-SDC4 were only observed between mouse germ cells and epithelial cells or granulosa cells.

We next reversed the signaling direction by designating germ cells as the signal receptor to examine incoming signals from surrounding somatic cell types (Figure 5C). In contrast to the analysis in which germ cells were treated as signaling sources, overall interaction strength was reduced, and only species-specific signaling pathways were identified. For example, WNT5A-FZD3 and KITL-KIT interactions, which are known to trigger PI3K/Akt signaling and promote fibroblast adhesion⁹⁰, and to play crucial roles in primordial follicle formation, activation, and follicular development⁹¹, respectively, were observed between bovine granulosa cells and germ cells. THBS1-CD47 signaling was detected specifically between mouse granulosa cells and germ cells. Interestingly, we did not find significant enrichment of ligand-receptor interaction pairs between human somatic cells and germ cells.

Cross-species cell type classification using a support vector machine (SVM) model

To identify conserved cell types across three species, we trained a support vector machine (SVM) classifier on the human dataset to distinguish major gonadal cell types and applied the trained model to predict mouse and bovine cell types (Figure 6). In mice, erythroid cells, immune cells, endothelial cells, neural cells, PV, and germ cells exhibited high prediction probabilities, close to 1, indicating strong cross-species conservation (Figure 6A). Mesenchymal, mesothelial, and epithelial cells showed moderately high probabilities (∼0.6-0.8), suggesting partial conservation. In contrast, granulosa cells showed low prediction probabilities (<0.2), suggesting pronounced cross-species divergence on these cell types (Figure 6A). In bovine, prediction patterns differed from those of the mouse (Figure 6B). For example, glomerular mesonephros cells, endothelial cells, immune cells, neural cells, PV, erythroid cells, and germ cells showed high probabilities (>0.8). Mesothelial, mesenchymal, epithelial, and granulosa cells showed low prediction probabilities (<0.4), with granulosa cells being the lowest (Figure 6B). These results indicate that some gonadal somatic cells such as immune cells exhibit more conserved transcriptomic profiles across species, as reflected by their higher cross-species prediction probabilities and greater overlap in UMAP space. In contrast, granulosa cells display clearer species-specific divergence, with distinct clustering observed across the three species (Figure 6C).

Figure 6.

To further characterize which genes contribute to granulosa cell divergence across species, we trained a second SVM using granulosa cells from all three species and extracted genes with the highest model weights (Figure 6D). In bovine granulosa cells, top-weighted genes, such as LAMA2⁹² and MMP14⁹³, were associated with extracellular matrix organization and remodeling of the cellular microenvironment. In human granulosa cells, genes such as RHOB⁹⁴ and ARHGAP28⁹⁵ were enriched, implicating Rho GTPase-mediated cytoskeletal signaling. In mouse granulosa cells, genes including SOX4⁹⁶ and FST⁹⁷ were prominent, which were related to TGF-β/Activin-related signaling. Together, these features highlight species-specific regulatory features in granulosa cells.

Finally, to investigate the identity of the unclassified bovine-specific cell population (Figure 1B), we trained an SVM on bovine data excluding this population and used the model to predict its closest cell type identity (Supplementary Figure 5A). The majority of cells were classified as mesenchymal or granulosa cells, together accounting for >90% of predictions, with smaller contributions from germ cells and epithelial cells. Compared with mesenchymal cells, the unclassified bovine-specific population showed a distinct pattern, such as increased expression of FXYD3⁹⁸, involved in Na⁺/K⁺-ATPase–mediated ion regulation, and HS3ST5⁹⁹, which functions in heparan sulfate sulfotransferase (Supplementary Figure 5B). Indeed, this population did not express canonical granulosa cell markers but instead partially expressed the mesenchymal lineage-associated marker PDGFRA¹⁰⁰ (Supplementary Figure 5C). As early ovarian somatic cells include precursors of multiple lineages, including steroidogenic cells, we next assessed the expression of steroidogenesis-related genes. These cells showed clear expression of CYP17A1¹⁰¹, and partial expression of genes involved in steroid production and transport, including INSL3¹⁰², HSB3D1¹⁰³, and STAR¹⁰⁴ (Supplementary Figure 5C). Together, these results suggest that the unclassified bovine-specific somatic population may exhibit molecular features consistent with a steroidogenic lineage.

Discussion

We constructed an integrated single-cell transcriptomics atlas of early female gonadal development across bovine, human, and mouse to systematically identify conserved and divergent features of gonadal differentiation and define underlying molecular regulatory mechanisms. Overall, we identified 11 shared gonadal cell types, as well as a bovine-specific cell population that may exhibit steroidogenic features. Epigenetic regulators exhibited dynamic changes in both germ cells and granulosa cells. Developmental trajectory analysis revealed common transcriptional dynamic pattern regulating germ cell meiotic entry and granulosa cell differentiation, alongside divergence in the top functional gene pathways and regulons in each species. Cell-cell communication analysis further demonstrated that early gonad development is supported by extensive and dynamic signaling interactions between germ cells and surrounding somatic cells. Finally, implementation of a support vector machine (SVM)-based framework enabled quantitative assessment of transcriptional conservation and provided a predictive model for cross-species cell lineage classification.

Classification of a bovine-specific cell type

In our analysis, we used an SVM-based classifier to identify a unique bovine-specific cell type, which characterized by steroidogenic features, including expression of CYP17A1 and HSD3B1. In humans, CYP17A1 expression was not detected in early female gonadal cell types between PCW5-12, however, it was observed in male Leydig and interstitial cells in the same time window¹⁰⁵. CYP17A1 is also a known marker of theca cells¹⁰⁶, which typically differentiate in the ovarian follicle around puberty in females. In mice, theca cells progenitors arise from two populations, including an ovary-intrinsic precursor population detectable as early as E10.5 in the genital ridge, and a mesenchymal population originating from the mesonephros that migrates into the developing ovary around E12.5¹⁰⁷. However, neither of these theca precursor populations exhibits strong steroidogenic gene expression (e.g., HSD3B1) during fetal development¹⁰⁷. In contrast, we identified a bovine-specific cell population expressing steroidogenic genes (CYP17A1 and HSD3B1) that presents from E38-E112, and does not align with the developmental patterns described in humans and mice. This divergence in both gene expression and developmental timing suggests that this population may represent either a bovine-specific early steroidogenic lineage in the female gonad or a theca precursor population with species-specific temporal state.

Conserved and divergent regulatory programs in germ cell and granulosa cell development

Our developmental trajectory analysis shows a conserved gene expression pattern during the transition from PGCs to meiotic entry across cattle, humans, and mice. This process begins with the downregulation of pluripotency factors POU5F1⁴⁶, NANOG⁴⁷, and SOX17⁴⁹, followed by the activation of maturation markers like KIT⁵⁰ and ELAVL2⁵¹. While SOX17 is known as a key specifier in humans⁴⁹, its consistent dynamic in our bovine data confirms its potential role for ungulate germline development as well. All three species eventually activate the same meiotic program, including structural genes SYCP1⁵², SYCP2⁵³, SYCP3⁵⁴, and the epigenetic regulator ZCWPW1⁵⁶. The shared timing of these genes across species suggests that the core mechanism for entering meiosis is highly stable in mammals.

Similarly, granulosa cell trajectories highlight a conserved role for the AP-1 family, specifically FOS and JUNB⁸³. These factors act as common regulators for follicle maturation across all three species. While a previous study showed they coordinate metabolic shifts and cholesterol synthesis in adult periovulatory granulosa cells¹⁰⁸, our results suggest they may also play a role in embryonic granulosa cells. Interestingly, PBX1 showed the strongest regulon activity in bovine germ cells and was also elevated in bovine granulosa cells compared to humans and mice. Since PBX1 is known to regulate genes involved in steroid production⁶³, its high activity in both cell types may be linked to the bovine-specific cell population we identified that expresses steroidogenic markers. This suggests that in cattle, PBX1 might help coordinate a specific endocrine environment during early gonadal development that is not as prominent in the other two species. FOXL2 also showed species divergence. As a critical granulosa cell determinant in humans and mice²¹, it was detected at late pseudotime in these species but not in bovine, suggesting a delayed granulosa cell maturation program in cattle, in contrast to its strong expression in adult bovine follicular granulosa cells²⁴.

Cell-cell communication between germ cells and gonadal somatic cells

Cell-cell signaling interaction analysis further provides insight into how the germline niche may be established during early gonadal development. Compared with the extensive communication among somatic cell populations, germ cells-associated signaling was relatively limited, with germ cells primarily acting as signaling sources rather than receivers. This pattern raises the possibility that early germ cells may contribute to shaping their local microenvironment rather than functioning solely as passive recipients of niche signals. Conserved ligand-receptor pairs were observed, including MDK–SDC4 in all species and NECTIN3–NECTIN2 in cattle and mice, suggesting that growth factor signaling and adhesion-mediated interactions may represent conserved mechanisms contributing to germ cell niche establishment. MDK shows prominent expression in human granulosa cells, and MDK-syndecan interactions have been reported to regulate primordial follicle formation in humans and mice¹⁰⁹. Similarly, nectin family adhesion molecules mediate cell–cell junctions in ovarian follicles¹¹⁰, and NECTIN2 has been observed to be expressed at cell-cell adhesion sites within the granulosa cell layer in the primary and preantral follicles in mice¹¹¹, suggesting that NECTIN3-NECTIN2 interactions may play a conserved role in mediating cell-cell adhesion during early gonadal development in cattle and mice.

Cross-species cell type classification and prediction

In addition, our cross-species cell type classification using an SVM revealed differential levels of transcriptional conservation among cell types. germ cells exhibit greater cross-species similarity, whereas granulosa cells display more pronounced cross-species divergence. This contrast likely reflects the fundamental and evolutionarily conserved roles of germ cells in early gonadal development¹¹². In contrast, granulosa cells’ functions are more closely associated with species-specific follicular development and reproductive physiology¹¹³.

The top-weighted genes identified by the SVM classifier for species-specific granulosa cell discrimination converge on genes related to extracellular matrix (ECM) organization and microenvironmental regulation. In cattle, highly ranked genes such as LAMA2⁹² and MMP14⁹³ are directly involved in ECM structure and remodeling. In humans, genes associated with Rho GTPase-mediated cytoskeletal signaling were enriched, a pathway closely linked to cell–matrix interactions and mechanical sensing¹¹⁴. In mice, the top genes are components of TGF-β/Activin signaling, which plays a key role in regulating ECM production and tissue remodeling¹¹⁵. Together, these results suggest that interspecies differences in granulosa cells are likely associated with divergence in ECM-related regulatory programs. Beyond revealing patterns of transcriptional conservation across species, SVM classifier also provides a predictive framework for cross-species cell type identification¹¹⁶. For example, by leveraging conserved transcriptional signatures, SVM classifier perform the best among 22 single cell classifiers in annotating cell populations in newly generated single-cell datasets, particularly in cases where canonical marker genes are incomplete, poorly characterized, or absent¹¹⁶. Such a strategy may be especially valuable for studies in non-model organisms, such as ruminant livestock, where cell type annotation often remains challenging¹¹⁷. Moreover, SVM classifier offers a scalable framework for integrating comparative single-cell datasets and for systematically identifying conserved and divergent cellular programs across species.

Despite these insights, several limitations of cross-species transcriptomic analysis should be considered. First, cross-species comparisons rely on approximate developmental stage alignment, and differences in sampling windows and dataset composition may influence the detection of conserved or species-specific cell populations. Second, our analysis is restricted to one-to-one orthologous genes, which limits the ability to capture species-specific genes or gene family expansions that may play important roles in lineage diversification. Future studies could address these limitations by expanding gene mapping strategies by incorporating all ortholog types and by integrating additional layers of information, such as regulatory or epigenetic features, to improve the robustness of cross-species comparisons.

Conclusion

In summary, our integrative cross-species single-cell analysis provides a comprehensive framework for understanding early ovarian development across mammals. We identified both conserved gene regulatory programs and species-specific features that shape female germline and gonadal somatic cell differentiation. Our findings reveal a high degree of conservation in the transcriptional programs underlying germ cell development and meiotic entry, supporting the conservation of core developmental mechanisms. In contrast, granulosa cells exhibit species-specific divergent extracellular matrix-related gene regulatory networks, indicating how evolutionary differences may shape reproductive strategies and developmental timing. Beyond these biological insights, our study establishes a generalizable framework for cross-species single-cell transcriptomic integration, quantitative assessment of transcriptional conservation, and cell type annotation. This resource lays the groundwork for future comparative studies in reproductive biology and offers a foundation for exploring cell-type-specific regulatory programs across species.

Acknowledgements

The authors acknowledge all members from the Tang lab and Wei lab for assistance with sample collection, and all members from the Duan lab for support with study design and data analysis.

Data availability

Human (PCW6, PCW7, PCW10, PCW12, PCW13, and PCW16) and mouse (E11.5, E12.5, E14.5, E16.5, and E18.5) single-cell RNA-seq datasets used in this study are publicly accessible from ArrayExpress (E-MTAB-10551) and GEO (GSE136441). Single-cell RNA sequencing data for fetal bovine gonads are available in the China National Center for Bioinformation (CNCB) database under accession number: OMIX012725.

Additional information

Funding

This study was funded by the National Science Foundation through Grant # 2213824 to J.E.D. and the IISAGE Consortium. This project was supported by Biological Breeding-National Science and Technology Major Project (2023ZD0407504) to Y.T., and Chinese Universities Scientific Fund (245-2023-F2010123001, 245-2025-Z1090224009) to Y.T.

Authors’ contributions

Y.F. performed data processing, integration, computational analyses, and figure preparation. S.H. collected the samples and conducted the experiments. J.Z., S.X., and Y.S. assisted in sample collection. J.E.D., T.Y., and Y.W. supervised the study and provided conceptual guidance. Y.F. drafted the manuscript, and all authors reviewed, revised, and approved the final version.

Funding

National Science Foundation (NSF) (2213824)

• Jingyue Duan

MOE | Chinese Universities Scientific Fund (245-2023-F2010123001)

• Young Tang

MOE | Chinese Universities Scientific Fund (245-2025-Z1090224009)

• Young Tang

Additional files

Supplementary figures.

Supplmentary table.