Integrative single-cell transcriptomics of sheep ovarian development reveals dynamic transcriptional programs relevant to reproductive traits.

B.Z. developed the computational pipeline, performed data analysis, and contributed to the writing of the article. K.D. carried out formal data analysis and validation. S.L. organized and visualized the data and was responsible for data curation. G.Z. contributed to the sample collection. W.L. was involved in data visualization. W.L. participated in the formal data analysis. Y.C. contributed to sample collection. F.W. provided experimental resources. H.Z. contributed to writing, review, and editing. D.G. was responsible for conceptualization and lead article writing, review, and editing. Y.Z. led the conceptualization, formal analysis, writing, review, and editing, and was responsible for funding acquisition and project administration. All authors reviewed and approved the final article.

To investigate gene expression and cellular diversity during ovarian development, we profiled single-cell transcriptomes across five developmental time points representing key developmental stages in sheep: prenatal (E90), pre-puberty (M3), post-puberty (M6), and adulthood (Y2, Y4) ( Figure 1 A). We generated five scRNA-seq datasets, each corresponding to a developmental stage, pooling samples from three biological replicates per time point. Our analysis yielded 61,649 single-cell transcriptomes, with an average of 12,329 cells per sample ( Figure S1 A, Table S1 ). After applying quality control measures (see STAR Methods ) , we classified these cells into 22 distinct cellular clusters (CLUs) using the Seurat unsupervised clustering 35 ( Figures 1 B and S1 B). Based on the expression of well-defined marker genes, we identified nine major ovarian cell types ( Figures 1 B and 1C). Specifically, CLUs 0, 1, 2, 3, and 4 were classified as stromal cells ( n = 30,812), the largest population, based on the expression of markers including DCN , COL1A1 , and COL1A2 . 13 , 36 , 37 The next largest population was granulosa cells ( n = 12,537), identified in CLUs 5, 6, 7, and 8 by high expressions of markers of FST , FSHR , and CDH2 , 38 , 39 which played crucial roles in supporting oocyte growth, maturation, and ovulation. We then identified 530 oocytes, represented by CLUs 9 and 10, by classical markers of ZP3 , 28 DAZL , 40 and OOEP . 24 Immune cells ( n = 2,915) were assigned to CLUs 11, 12, and 13, based on the expression of PTPRC and CD53 , 20 , 38 while epithelial cells ( n = 2,380) were identified in CLUs 14, 15, and 16 by the expression of KRT19 and ALDH1A2 . 41 , 42 CLUs 17 and 18 represented endothelial cells ( n = 4,956), marked by high expression levels of CDH5 , CYYR1 , FLI1 , and PECAM1 . 24 , 43 , 44 Theca cells ( n = 4,013) were identified in CLU 19 by markers CYP11A1 and LHCGR , which were associated with follicle development and maturation. 45 The CLU 21, which exhibited the high expression of MYH11 , MYL9 , and TAGLN , 46 , 47 , 48 , 49 was identified as smooth muscle cells ( n = 3,226). Lastly, a small population of mesenchymal cells ( n = 280) was assigned to CLU 20, based on the expression of S100A4 . 50 Figure 1 Single-cell transcriptomic profiles across ovarian developmental stages (A) Illustration depicting the sampling of sheep ovaries at three developmental stages: prenatal (embryonic day 90 [E90)], youth (3 months [M3] and 6 months old [M6]), and adulthood (2 years [Y2] and 4 years old [Y4]). Corresponding human developmental stages are shown for comparison on the left side. (B) Uniform Manifold Approximation and Projection (UMAP) plot showing the distribution of 61,649 ovarian cells across 5 time points derived from 3 developmental stages (E90, M3, M6, Y2 and Y4) based on single-cell gene expression. The dimension reduction analysis classified ovarian cells into 21 clusters, indicated by different colors, which are further categorized into 9 cell types based on well-defined marker genes (see details in main text). (C) Violin plots depict the normalized gene expression of known cell-type gene markers that define 9 ovarian cell types. The color of each violin corresponds to a specific cell type, consistent throughout the article. (D) Proportion of 9 ovarian cell types (y axis) across 5 time points (E90, M3, M6, Y2 and Y4) (x axis). (E) Proportion of cells in predicted cell cycle phases across 9 ovarian cell types in 3 developmental stages. G1 (Gap 1 Phase): The initial phase post-cell division where the cell grows, synthesizes mRNA and proteins essential for DNA synthesis, and prepares for DNA replication; S (Synthesis Phase): DNA replication occurs, ensuring each new cell receives an identical set of DNA; G2M (Gap 2 to Mitosis Phase): The final preparation stage where cells grow and ready themselves for mitosis, completing any needed repairs and synthesis before division. Single-cell transcriptomic profiles across ovarian developmental stages (A) Illustration depicting the sampling of sheep ovaries at three developmental stages: prenatal (embryonic day 90 [E90)], youth (3 months [M3] and 6 months old [M6]), and adulthood (2 years [Y2] and 4 years old [Y4]). Corresponding human developmental stages are shown for comparison on the left side. (B) Uniform Manifold Approximation and Projection (UMAP) plot showing the distribution of 61,649 ovarian cells across 5 time points derived from 3 developmental stages (E90, M3, M6, Y2 and Y4) based on single-cell gene expression. The dimension reduction analysis classified ovarian cells into 21 clusters, indicated by different colors, which are further categorized into 9 cell types based on well-defined marker genes (see details in main text). (C) Violin plots depict the normalized gene expression of known cell-type gene markers that define 9 ovarian cell types. The color of each violin corresponds to a specific cell type, consistent throughout the article. (D) Proportion of 9 ovarian cell types (y axis) across 5 time points (E90, M3, M6, Y2 and Y4) (x axis). (E) Proportion of cells in predicted cell cycle phases across 9 ovarian cell types in 3 developmental stages. G1 (Gap 1 Phase): The initial phase post-cell division where the cell grows, synthesizes mRNA and proteins essential for DNA synthesis, and prepares for DNA replication; S (Synthesis Phase): DNA replication occurs, ensuring each new cell receives an identical set of DNA; G2M (Gap 2 to Mitosis Phase): The final preparation stage where cells grow and ready themselves for mitosis, completing any needed repairs and synthesis before division. Although the same cell types were present at all developmental stages, their proportions varied dynamically ( Figures 1 D, S1 C and S1D). Notably, granulosa and stromal cell populations decreased rapidly after M6, while endothelial, smooth muscle, and immune cells increased ( Figure 1 D). The further cell cycle analysis indicated that most M3 and M6 cells were enriched in the G2/M phase, suggesting active cell proliferation and preparation for mitosis, which aligned with the critical stages of follicular development and ovarian growth during pre-puberty and post-puberty. Conversely, oocytes remained quiescent state (G1) at the stages Y2 and Y4, reflecting post-pubertal follicular reserve maintenance ( Figure 1 E). To study gene expression across cell types, we detected cell type-specifically expressed genes by comparing a target cell type with the rest ( Table S2 ). As shown in Figure 2 A, for instance, the FST gene, a key folliculogenesis regulator, 51 was highly expressed in granulosa cells. In the developing mouse ovary, its expression was regulated by FOXL2 (forkhead-domain transcription factor L2), 52 which was also highly expressed in granulosa cells in this study ( Table S2 ). As expected, these cell type-specifically expressed genes reflect the cell-type identity and physiology. Similarly, genes such as FSHR and CDH2 were specifically and highly expressed in granulosa cells, whereas CD44 , CD52 , and CD74 were mainly expressed in immune cells. Muscle-related genes, including MYH11 and MYL6, were predominant in smooth muscle cells. The expression of FSHR , a receptor essential for follicle growth and ovulation, further supports the physiological role of granulosa cells. 53 The functional relevance of these cell-type-specific genes was supported by functional enrichment analysis, which showed significant associations (FDR <0.05) with known biological pathways ( Figures 2 B and S2 ; Table S3 ). For example, genes specific to endothelial cells were significantly enriched in pathways related to blood vessel morphogenesis and endothelium development. Similarly, immune cell-specific genes were associated with immune response and lymphocyte activation, while muscle-specific genes were enriched for muscle development and contraction. To further elucidate transcriptional regulation within cell types, we employed the SCENIC tool 54 and identified key transcription factors (TFs) associated with each cell type. In granulosa cells, for example, FOXO1, HMGA2, PPARG, ESR1, and GATA4 were identified. Notably, more than 90% of their target genes were cell-type-specific ( Figure 2 C), indicating that these TFs drive distinct ovarian cell identities and functions. Figure 2 Identification of cell-type-specifically expressed genes (A) Heatmap shows the relative expression of specifically expressed genes in each cell type. Expression level was shown with its Z score. (B) Gene Ontology (GO) terms significantly enriched (False Discovery Rate, FDR <0.05) by cell type-specific genes. (C) TF-target gene pairs in granulosa cells. The light green dots show the granulosa cell specifically expressed genes, while light purple dots show non-cell-type specifically expressed genes. (D) Heatmap showing the scaled expression patterns of hotspot genes for POF/POI 55 , 56 , 57 (right) in different cell types. POF: premature ovarian failure; POI: primary ovarian insufficiency. Identification of cell-type-specifically expressed genes (A) Heatmap shows the relative expression of specifically expressed genes in each cell type. Expression level was shown with its Z score. (B) Gene Ontology (GO) terms significantly enriched (False Discovery Rate, FDR <0.05) by cell type-specific genes. (C) TF-target gene pairs in granulosa cells. The light green dots show the granulosa cell specifically expressed genes, while light purple dots show non-cell-type specifically expressed genes. (D) Heatmap showing the scaled expression patterns of hotspot genes for POF/POI 55 , 56 , 57 (right) in different cell types. POF: premature ovarian failure; POI: primary ovarian insufficiency. To determine whether genes associated with ovarian diseases exhibit cell-type-specific expression, we analyzed genes linked to premature ovarian failure (POF) and primary ovarian insufficiency (POI). 55 , 56 , 57 Our findings revealed distinct expression profiles across different ovarian cell types ( Figure 2 D). For instance, FST , FSHR , ESR1 , and CDH2 were highly expressed in granulosa cells, while GDF9 , POU5F1 , and NANOS3 were predominantly expressed in oocytes. The GDF9 , encoding growth differentiation factor-9, plays an essential role in folliculogenesis, oogenesis, and ovulation. 58 Interestingly, genetic polymorphisms in the GDF9 gene have been strongly associated with fecundity in various sheep breeds. 59 , 60 , 61 These findings underscore the importance of specific ovarian cell types in maintaining ovarian homeostasis, with substantial implications for fertility and reproductive performance. The inclusion of multiple developmental stages in our study enabled us to explore the temporal dynamics of gene expression across ovarian cell types. To investigate how intercellular communication evolves over time, we utilized CellChat 62 to analyze cell-cell interactions between different cell types ( n = 9) across developmental stages. Our analysis revealed a dynamic pattern of intercellular signaling, with notable variations in the number and strength of interactions throughout ovarian development ( Figure 3 A). One key dynamically regulated pathway was the Notch signaling pathway, which plays critical roles in regulating cell proliferation, differentiation, and apoptosis, and follicular development 63 ( Figure S3 ). Notably, in adult ovaries, oocytes demonstrated minimal detectable interactions with other cell types, possibly reflecting a reduced intercellular communication associated with aging. 19 , 20 In contrast, at young stages (i.e., M3 and M6), granulosa and theca cells exhibited heightened intercellular interactions, both within and between their respective cell types ( Figure 3 B). This elevated communication aligns with their critical roles in supporting follicular growth and maturation during puberty. 64 Figure 3 Dynamic gene expression of sheep ovarian cells (A) The number of interactions for each cell type in the developmental stages. (B) Circle plot of the ligand–receptor pairs of all cell types, with the thickness of the string representing the weight of the ligand–receptor pairs. (C) Z score scaled expression levels of developmental stage-specific genes (y axis) across 5 time points (x axis) (left), with corresponding enriched Gene Ontology (GO) terms (middle) and transcription factor (TF) binding motifs (right). Dynamic gene expression of sheep ovarian cells (A) The number of interactions for each cell type in the developmental stages. (B) Circle plot of the ligand–receptor pairs of all cell types, with the thickness of the string representing the weight of the ligand–receptor pairs. (C) Z score scaled expression levels of developmental stage-specific genes (y axis) across 5 time points (x axis) (left), with corresponding enriched Gene Ontology (GO) terms (middle) and transcription factor (TF) binding motifs (right). To further investigate stage-specific expression patterns, we identified genes uniquely expressed at each developmental stage ( Table S4 ). In granulosa cells, for example, we identified 1,821 genes uniquely expressed prior to birth ( Figures 3 C; Table S4 ), many of which were enriched in metabolic and hormone signaling pathways. Among them, JUND, a key component of AP1-DNA binding complex, was the most highly enriched transcription factor essential for the maturation and differentiation of granulosa cells. 65 In postnatal stages (M3 and M6), these stages featured distinct gene sets associated with hormone secretion and ovarian follicle migration, reflecting the ovary’s transition toward follicle maturation and ovulation. Notable TFs include FOXO1, which in sheep granulosa cells is involved in proliferation inhibition, promoting apoptosis, cell cycle progression, and steroidogenesis. 66 Another key TF, ESR1 (estrogen receptor 1), played a crucial role in dominant follicle development in ruminants, distinct from its function in polyovulatory species such as mice. 67 At adult stages (Y2 and Y4), FOXL2, a member of the forkhead box family, was the most enriched TF associated with folliculogenesis, estrogen production, granulosa cell proliferation, and tumorigenesis. 68 Another key TF, EGR1 (early growth response 1), a zinc finger TF regulates apoptosis through fibroblast growth factor signaling in granulosa cells. 69 The expression of EGR1 in aged granulosa cells suggests an enhancement of apoptotic pathways, potentially contributing to age-associated ovarian decline, which correlates with the observed loss of oocyte communications. Given the pivotal role of granulosa cells in ovarian function, particularly during young reproductive stages, we conducted unsupervised clustering analysis on granulosa cells ( n = 12,537), identifying nine granulosa subclusters. These subclusters were categorized into preantral, antral, and atretic granulosa cell subtypes based on marker gene expression 20 ( Figures 4 A and 4B; Table S5 ). Preantral granulosa subtypes (subcluster 0, 1, and 2) showed high expressions of IGFBP5 , 25 GATM , 13 and COL18A1 20 ; antral granulosa subtypes (subcluster 3, 4, 5, and 6) were marked by INHBB , 70 FST 33 and GJA1 20 ; and atretic granulosa subtypes (subcluster 7 and 8) by ITIH5 and GHR 20 ( Figure 4 B). The composition of these granulosa cell subtypes varied dynamically across developmental stages ( Figure 4 C). A higher proportion of preantral granulosa cells was observed in the embryonic stage (E90) and adult stages (Y2, Y4), likely due to reduced follicle recruitment and slower follicular turnover. In contrast, the young stages (M3, M6) exhibited a lower proportion of preantral granulosa cells, reflecting active recruitment of follicles toward the antral stage, which depleted the preantral granulosa pool. Accordingly, we observed an exceptionally higher proportion of antral granulosa cell subtypes in the young stages (M3, M6), and the population of atretic granulosa cell subtype was exclusively present in these stages. These observations might reflect high follicular turnover where many follicles initiate growth, but only a few reach ovulation. 71 Figure 4 Characterization of granulosa cell type diversity (A) Uniform Manifold Approximation and Projection (UMAP) analysis of granulosa cells (GCs), whose subclusters were grouped into preantral, antral, and atretic granulosa cell subtypes based on marker gene expressions (see details in the main text). (B) Heatmap depicting marker gene expression (z-scored) used for granulosa cell subtype identification. (C) The proportion of granulosa cell subtypes at 5 time points across 3 developmental stages. E90: embryonic 90 days; M3: 3 months old; M6: 6 months old; Y2: 2 years old; Y4: 4 years old. Color legend is the same as the panel (B). (D) Monocle2 pseudotime analysis. Scatterplots showing the cell trajectories by modeled pseudotime (left), granulosa cell subtypes (middle), predicted cell states (right). (E) Heatmap showing pseudotime gene clusters differentially expressed during granulosa cell fate commitment (left), the top significant ( p < 0.05) GO terms and KEGG pathways for each gene set (middle), and expression levels of representative genes in cells ordered along the pseudotime trajectory (colored based on states to which they are assigned in the panel (D), right). Characterization of granulosa cell type diversity (A) Uniform Manifold Approximation and Projection (UMAP) analysis of granulosa cells (GCs), whose subclusters were grouped into preantral, antral, and atretic granulosa cell subtypes based on marker gene expressions (see details in the main text). (B) Heatmap depicting marker gene expression (z-scored) used for granulosa cell subtype identification. (C) The proportion of granulosa cell subtypes at 5 time points across 3 developmental stages. E90: embryonic 90 days; M3: 3 months old; M6: 6 months old; Y2: 2 years old; Y4: 4 years old. Color legend is the same as the panel (B). (D) Monocle2 pseudotime analysis. Scatterplots showing the cell trajectories by modeled pseudotime (left), granulosa cell subtypes (middle), predicted cell states (right). (E) Heatmap showing pseudotime gene clusters differentially expressed during granulosa cell fate commitment (left), the top significant ( p < 0.05) GO terms and KEGG pathways for each gene set (middle), and expression levels of representative genes in cells ordered along the pseudotime trajectory (colored based on states to which they are assigned in the panel (D), right). To model the developmental trajectory of granulosa cells, we utilized Monocle2 software 72 to conduct pseudotime analyses, which revealed five unique cell states transitioning through two major branch points ( Figure 4 D). In alignment with the shifting proportions of granulosa cell subtypes over time ( Figure 4 C), a higher proportion of preantral granulosa cells transitioned into antral subtypes through the branchpoint 1, whereas at later stages, more cells progressed toward atretic subtypes via the branchpoint 2 ( Figure 4 D). Specifically, cell state 3 transitioned into states 4 and 2 through the branchpoint 1, while state 2 further differentiated into states 5 and 1 through the branchpoint 2. At the branchpoint 1, where preantral and antral granulosa cells predominated, we identified three major gene clusters ( Figure 4 E). Cluster 1: Early-stage genes ( n = 541), enriched in GTPase activity and AMPK/FoxO signaling pathways, which are essential for primordial follicle activation. Cluster 2: Genes highly expressed during young stages ( n = 354), associated with hormone secretion (e.g., FST, FSHR, INHBA, HIF1A, PPARG, EGFR ). Cluster 3: Adult-stage genes ( n = 1,032), enriched in protein targeting and metabolic regulation, reflecting the shift from follicular growth to atresia. A similar pipeline was applied to immune cells, revealing three major immune subtypes, including T lymphocytes, macrophages, and neutrophils ( Figure S5 ). The macrophages were identified by the expression of CD74 , AIF1 and C1QA , 28 , 73 the neutrophils were enriched for CSF3R , S100A8 and S100A9 , 74 whereas the T lymphocytes showed high expression of CD3G 75 ( Figure S5 A). The proportion of macrophages decreased while T lymphocytes increased from embryonic to adulthood stages ( Figure S5 B). CellChat analysis further showed strong interactions between macrophages and granulosa cells, suggesting a key role for immune cells in follicular dynamics ( Figure S5 C). Furthermore, a small subset of Natural Killer T (NKT) cells ( n = 61) was identified using the canonical marker gene CD3D 8 ( Figures S5 D and S5E). To explore cross-species similarities, we integrated human ovary scRNA-seq data 76 , 77 ( Figure S6 , Table S6 ) with our sheep dataset, aligning 80,912 human ovarian cells with 61,649 sheep ovarian cells based on 13,639 one-to-one orthologous genes ( Figure 5 A). While there was substantial overlap in ovarian cell types between sheep and humans, our integrative cross-species analysis identified an additional mesenchymal cell population in sheep that was not detected in single-species analysis ( Figure 5 A). Mesenchymal cells are known to promote ovarian function recovery by inhibiting granulosa cell apoptosis and follicular atresia, potentially through the upregulation of anti-Müllerian hormone and follicle-stimulating hormone receptor expression in granulosa cells. 78 Further analysis revealed conserved cell population frequencies across comparable developmental stages between sheep and humans ( Figure S6 C, Figure 1 A). Specifically, a higher proportion of granulosa cells was detected during embryonic and early postnatal stages, while smooth muscle cells increased with aging, reflecting age-related ovarian remodeling. However, stromal cells displayed divergent patterns: they were more abundant in embryonic sheep ovaries, whereas in humans, their abundance increased postnatally. Furthermore, we investigated gene expression conservation by conducting Meta-neighbor analysis using the MetaNeighbor tool. 79 Overall, we found that gene expression was conserved between sheep and humans, with an Area Under the Receiver Operating Characteristic (AUROC) score greater than 0.75 in a majority of cell type pairs ( Figure 5 B). Functional enrichment of orthologous genes confirmed biological relevance, with conserved cell types showing expected physiological functions. For example, macrophage-enriched genes were associated with leukocyte degranulation, neutrophil mediated immunity ( Figure 5 C). Figure 5 Comparative analysis of single-cell transcriptomes from human and sheep ovaries (A) Uniform Manifold Approximation and Projection (UMAP) showing the locations of cell types of ovine ovaries generated in this study and human ovaries from Chitiashvili et al. 76 and Wu et al. 77 (B) The area under the receiver operating characteristic curve (AUROC) calculated using 1-to-1 orthologous genes ( n = 13,639) of the same cell types in sheep and humans. The asterisks indicated the AUROC >0.75, a threshold used in this study to justify the conservation of gene expression in the same cell type. Human (rows) and sheep (columns) cell types were separately clustered based on their AUROC values, using the Euclidean method 80 to compute distances between cell types. The hierarchical tree was constructed using the complete linkage method. (C) Functional enrichment analysis of orthologous genes in conserved cell types between human and sheep. Regulon specificity score (RSS) of transcription factors (TFs) in macrophages (D), T lymphocytes (E), granulosa cell (F). (G) Cell-cell communication networks in the ovaries of sheep (left) and humans (right). The line width indicates the strength of cell-cell interactions, with thicker lines representing stronger interactions. Comparative analysis of single-cell transcriptomes from human and sheep ovaries (A) Uniform Manifold Approximation and Projection (UMAP) showing the locations of cell types of ovine ovaries generated in this study and human ovaries from Chitiashvili et al. 76 and Wu et al. 77 (B) The area under the receiver operating characteristic curve (AUROC) calculated using 1-to-1 orthologous genes ( n = 13,639) of the same cell types in sheep and humans. The asterisks indicated the AUROC >0.75, a threshold used in this study to justify the conservation of gene expression in the same cell type. Human (rows) and sheep (columns) cell types were separately clustered based on their AUROC values, using the Euclidean method 80 to compute distances between cell types. The hierarchical tree was constructed using the complete linkage method. (C) Functional enrichment analysis of orthologous genes in conserved cell types between human and sheep. Regulon specificity score (RSS) of transcription factors (TFs) in macrophages (D), T lymphocytes (E), granulosa cell (F). (G) Cell-cell communication networks in the ovaries of sheep (left) and humans (right). The line width indicates the strength of cell-cell interactions, with thicker lines representing stronger interactions. We next examined conservation in ovarian gene regulatory networks by comparing transcription factor regulon activity between the two species. Regulon specificity score (RSS) of TFs showed significant correlations across several key ovarian cell types, including: ( Figures 5 D–5F and S7 ) macrophages ( r = 0.57, FDR = 3.35 × 10 −16 , Figure 5 F), T lymphocytes ( r = 0.38, FDR = 2.36 × 10 −7 , Figure 5 D), and granulosa cell ( r = 0.31, FDR = 6.42 × 10 −5 , Figure 5 F). We also observed notable similarities in cellular communication and signaling pathways between the two species, as shown in Figure 5 G. These analyses collectively indicated there were considerable conservation in cell types, cell type frequency, gene expression pattern and networks between human and sheep developing ovaries. To further explore whether the cross-species transcriptional conservation could provide insights into the genetic basis of reproductive traits and diseases, we analyzed genome-wide association studies (GWASs) data for nine human complex reproductive traits and three diseases ( Table S7 ). Polygenic regression analysis 81 revealed significant associations between sheep ovarian cell types and human reproductive traits. Notably, T lymphocyte (FDR = 2.48 × 10 −10 ) and macrophage (FDR = 2.15 × 10 −7 ) were strongly associated with fetal post-term birth, and granulosa cells showed significant associations with fetal pre-term birth (FDR = 4.27 × 10 −8 ) and endometriosis (FDR = 6.24 × 10 −8 ) ( Figure 6 A). Focusing on PCOS as a case study, trait-relevant score (TRS) 81 analysis identified smooth muscle cells as the most relevant cell type, followed by T lymphocytes ( Figures 6 A and 6B). Further investigation of PCOS-associated TF highlighted the Jun family (e.g., JUN, JUND, JUNB) and FOS as key regulators ( Figure 6 C). These findings align with previous studies, 82 , 83 , 84 reinforcing their role in PCOS pathophysiology. Figure 6 Integration of sheep ovary single-cell transcriptomes with human reproductive traits and disease (A) Associations (i.e., -log 10 FDR) of sheep ovarian cell types with human reproductive traits and diseases by using scPagwas. 81 (B) Trait-relevant score (TRS, left) and the corresponding significance (-log 10 p -value, right) for the correlations of conserved sheep ovarian cell types with polycystic ovary syndrome (PCOS). The red line represents the significant level. (C) PCOS-relevant genes ranked by the Pearson correlation coefficient (PCC) using scPagwas across all individual cells. (D) Associations (i.e., -log 10 FDR) of smooth muscle cell T lymphocyte with PCOS across 5 time points derived from 3 developmental stages (E90, M3, M6, Y2 and Y4) by using scPagwas. 81 (E) KEGG pathways relevant to PCOS in smooth muscle cells and T lymphocytes across developmental stages. KEGG: Kyoto Encyclopedia of Genes and Genomes. Integration of sheep ovary single-cell transcriptomes with human reproductive traits and disease (A) Associations (i.e., -log 10 FDR) of sheep ovarian cell types with human reproductive traits and diseases by using scPagwas. 81 (B) Trait-relevant score (TRS, left) and the corresponding significance (-log 10 p -value, right) for the correlations of conserved sheep ovarian cell types with polycystic ovary syndrome (PCOS). The red line represents the significant level. (C) PCOS-relevant genes ranked by the Pearson correlation coefficient (PCC) using scPagwas across all individual cells. (D) Associations (i.e., -log 10 FDR) of smooth muscle cell T lymphocyte with PCOS across 5 time points derived from 3 developmental stages (E90, M3, M6, Y2 and Y4) by using scPagwas. 81 (E) KEGG pathways relevant to PCOS in smooth muscle cells and T lymphocytes across developmental stages. KEGG: Kyoto Encyclopedia of Genes and Genomes. Focusing on different developmental stages, we observed that T lymphocytes showed the strongest association with PCOS during adult stages (Y2 and Y4), consistent with previous reports 21 ( Figure 6 D). In contrast, the relevance of smooth muscle cells observed in this study exhibited a significant but weaker association with PCOS, primarily during embryonic and adult stages ( Figure 6 D). At the embryonic stage, renin secretion was the only enriched pathway in smooth muscle cells associated with PCOS, whereas in the adult stage, pathways such as ubiquitin-mediated proteolysis and the AGE-RAGE signaling pathway in diabetic complications were significantly enriched, indicating potential mechanisms underlying PCOS progression ( Figure 6 E).

Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Yanli Zhang ( [email protected] ). This study did not generate new unique reagents. • Data that support the findings of this study are available in the key resources table and from the lead contact author upon reasonable request. • Codes for analyzing the single-cell data are available from the GitHub depository https://github.com/BingruZhao/Singlecell-RNA-seq-of-sheep-ovarian . • Any additional information required to reanalyze the data reported in this article is available from the lead contact upon request. Data that support the findings of this study are available in the key resources table and from the lead contact author upon reasonable request. Codes for analyzing the single-cell data are available from the GitHub depository https://github.com/BingruZhao/Singlecell-RNA-seq-of-sheep-ovarian . Any additional information required to reanalyze the data reported in this article is available from the lead contact upon request.

Through scRNA-seq analysis of sheep ovary tissue spanning prenatal to postnatal stages, we identified nine major ovarian cell types: stromal cells, granulosa cells, oocytes, immune cells, epithelial cells, endothelial cells, theca cells, mesenchymal cells, and smooth muscle cells. These cell types align with those consistently reported in previous studies across multiple species. 12 , 16 , 23 , 30 , 85 , 86 When comparing our findings to the previous single-cell study of the sheep ovary, 33 we observed substantial overlap in core cell types, including stromal cells, granulosa cells, immune cells, and endothelial cells. However, our study identified five additional cell types—oocytes, epithelial cells, theca cells, mesenchymal cells, and smooth muscle cells—whereas the previous study identified perivascular cells, a heterogeneous population of mesenchymal progenitors, underscoring the impact of cellular heterogeneity on cell-type annotation. The main difference may also stem from our sampling across multiple developmental stages, while the previous study focused on a single stage. 23 , 26 , 30 , 86 Notably, oocytes were exclusively detected at the E90 and Y2 samples, with their absence in M3, M6, and Y4 possibly reflecting sampling variability, technical artifacts, or physiological differences between developmental stages. Additionally, we identified six distinct subtypes within both the granulosa and immune cell populations, further illustrating the complexity of ovarian cell diversity. The immune cell population, in particular, exhibited significant complexity, consistent with prior pan-tissue analyses that identified over 100 immune cell subtypes, 87 suggesting that our study may still underestimate immune cell diversity, warranting further investigation. A particularly striking finding was the high proportion of immune cells at the Y4 stage, particularly T lymphocytes. This observation aligns with previous findings in mice, 19 , 20 suggesting that ovarian aging processes might initiate earlier than previously recognized, potentially beginning at the adult stage. The reduced oocyte-cell interactions and upregulation of apoptotic genes at Y4 further support this hypothesis. Although the same cell types were observed across developmental stages, their proportions shifted dynamically, reflecting the changing functional demands of the ovary. Granulosa and theca cells were most abundant at young stages (M3, M6), corresponding to an increased need for follicular growth and maturation during puberty. 88 In contrast, the adult stages (Y2, Y4) were characterized by a decline in follicle recruitment and reduced cell-cell communication, indicative of the onset of ovarian aging. 19 , 20 Furthermore, we compared gene expressions across developmental stages and explored how gene regulatory networks evolve to support distinct physiological processes. Our study identified stage-specific gene clusters linked to metabolic processes during the embryonic stage and hormone secretion and ovarian follicle migration during youth. Key TFs such as FOXO1, ESR1, FOXL2, and EGR1 appear to regulate stage-specific functions, including hormone production, cell proliferation, and apoptosis. 66 , 67 , 68 , 69 Interestingly, some TFs had both cell-type-specific and stage-specific functions, such as FOXO1 and ESR1, indicative of their critical roles in ovary function. For instance, positive associations were observed between ESR1 gene polymorphisms and PCOS across multiple populations. 89 Likewise, the FOXO1 was reported to be involved in the pathogenesis of PCOS. 90 These findings suggest that these TFs serve as potential biomarkers and therapeutic targets for reproductive disorders. 55 , 56 , 57 Our comparative analysis of human and sheep ovarian single-cell transcriptomes reveals strong conservation in cell types and gene expression patterns. A key finding was the identification of mesenchymal cells in sheep, likely due to enhanced resolution gained through dataset integration. 91 Using MetaNeighbor, 79 we observed significant similarity in cell-type pairs between humans and sheep, with an AUROC >0.75 across most cell type pairs ( Figure 5 B). However, our threshold was lower than AUROC >0.9 observed in a seven-species comparative study, including three vertebrates (human, mouse, and zebrafish) and four invertebrates ( C. intestinalis, C. elegans, Schmidtea, and Nematostella ). 92 This discrepancy may stem from the inclusion of diverse developmental stages, which introduced dynamic gene expression shifts ( Figure 3 C). These findings underscore the need to consider developmental context when conducting cross-species single-cell analyses. Indeed, the context-dependent nature of transcriptomic data substantially impacts the genetic architecture of complex traits. 93 For instance, our analysis indicated that T lymphocytes were predominantly relevant for PCOS at the adult stage—a finding consistent with previous studies 21 —we also observed that smooth muscle cells played a significant role in PCOS etiology during embryonic and adult stages, particularly through interactions with insulin metabolism-related pathways. These findings underscore the importance of profiling transcriptomes across developmental stages to gain a comprehensive understanding of the cellular physiology underlying complex reproductive traits and diseases, ultimately informing more effective therapeutic strategies. In summary, this study presents a comprehensive single-cell transcriptomic analysis of sheep ovaries, identifying nine major cell types and uncovering both cell-type-specific and developmental stage-specific gene expression patterns, potentially regulated by key TFs. Our findings defined dynamic shifts in ovarian cell composition and gene regulatory networks, identified key TFs involved in ovarian function and reproductive disorders, and demonstrated both conserved and species-specific transcriptional programs between sheep and humans, reinforcing the importance of developmental context in comparative analyses. By leveraging these insights, our study enhances the utility of sheep as a model organism for human reproductive research, providing a foundation for future studies on ovarian biology and related diseases. While this study advances our understanding of ovarian development in sheep and its relevance to human reproductive biology, several limitations should be acknowledged: 1) Limited species comparison: Our study focused on only two species (human and sheep), which may limit the generalizability of our findings. Expanding comparisons to additional mammalian models could provide broader evolutionary insights. 2) Lack of integration with other reproductive tissues: We analyzed only ovarian tissue, whereas reproductive processes involve multiple organs (e.g., uterus, hypothalamus). Future studies incorporating multi-organ transcriptomic data will provide a more holistic view of reproductive physiology. 3) Limited developmental time points: While we profiled five key stages, there remain gaps in our understanding of gene expression transitions across early embryogenesis and late reproductive aging. A more extensive temporal analysis would improve our resolution of dynamic ovarian changes. 4) Challenges in cross-species knowledge transfer: While our cross-species transcriptomic comparisons provide valuable insights, further refinement of computational methodologies is needed to improve alignment accuracy. Future efforts should integrate functional validation studies to confirm shared molecular mechanisms across species. Addressing these limitations in future research will be crucial for advancing reproductive biology and enhancing the translational potential of model organisms for studying human ovarian function and disease.

Sheep ( Ovis aries ), domesticated around 11,000 years ago based on archaeological and genetic evidence, 1 have adapted to diverse environments worldwide. 2 , 3 They play an essential role in agriculture, providing meat, wool, fur, and milk, with over 200 recognized breeds exhibiting significant variation in reproductive traits such as lambing frequency, litter size, and breeding season adaptability. Globally, sheep contribute approximately 8 million tons of meat annually ( https://www.fao.org/ ). Beyond their agricultural value, sheep serve as valuable biomedical models for studying various human disorders. 4 , 5 , 6 For example, due to similarities in lung size and structure, sheep are used in respiratory disease research. 5 In reproductive biology, sheep surpass rodents as models because their larger body size allows for precise ultrasound monitoring of ovarian follicular development, examination of placental development, and collection of hormonal and neurotransmitter profiles. 6 , 7 , 8 , 9 Despite differences in gestational length between humans (∼280 days) and sheep (∼150 days), key reproductive events such as gametogenesis, fertilization, implantation, parturition, and menstrual cycles are highly physiologically conserved. 7 Consequently, sheep have been widely utilized to study various human reproductive diseases and drug development. 6 For instance, polycystic ovary syndrome (PCOS), a common endocrine disorder characterized by hormonal imbalances, anovulation, and ovarian cysts formation, shares many physiological features between human and sheep ovaries. As such, sheep are effective models for studying PCOS pathophysiology and potential therapeutic inteventions. 7 The ovary is an essential organ for reproduction in mammals, comprising two main layers: the outer ovarian cortex and the inner ovarian medulla. 10 , 11 The cortex contains numerous follicles at various stages of development, each housing a primary oocyte surrounded by granulosa and theca cells. The inner medulla, consisting of loose connective tissue, blood vessels, lymphatic vessels, and nerves, supports the cortex but typically lacks follicles. 12 , 13 , 14 , 15 , 16 , 17 , 18 Within both the cortex and medulla, there are various types of immune cells, such as macrophages, dendritic cells, and T and B lymphocytes, contributing to tissue remodeling, follicular development, ovulation, and corpus luteum formation and regression. 19 , 20 Disruptions in ovarian structure or immune function can lead to disorders such as PCOS, ovarian cancer, ovarian cysts, and endometriosis. For example, T cell dysfunction has been implicated in PCOS-related infertility. 21 Thus, a comprehensive characterization of ovarian cell type and composition, along with gene transcription profiles, is essential for understanding reproductive mechanisms and developing targeted therapies for ovarian diseases. Single-cell RNA sequencing (scRNA-seq) has revolutionized our understanding of cellular heterogeneity and gene expression dynamics at single-cell resolution. 22 This technology has been extensively used to characterize ovarian cell types, compositions, gene transcription, and cell-cell communication across various species, including humans, 14 , 15 , 17 , 23 non-human primates, 18 , 24 mice, 12 , 18 , 25 , 26 rats, 27 yaks, 28 , 29 goats, 30 and Drosophila . 31 , 32 These studies have identified key ovarian cell populations, including oocytes, granulosa, endothelial, stromal, and immune cells. In sheep, Ge et al. 33 conducted a single-cell transcriptomic analysis of Hu sheep ovaries, a breed known for high fecundity, and highlighted the pivotal role of granulosa cells in prolificacy. However, no scRNA-seq studies have comprehensively examined ovarian development in sheep across multiple life stages, leaving the temporal dynamics of ovarian cell types and gene expression largely unexplored. Furthermore, most existing studies focus on single-species analyses, lacking cross-species comparisons that are essential for identifying conserved transcriptional programs and their implications for reproductive traits and disease. This knowledge gap limits the translational potential of sheep as a reproductive model, restricting its utility for studying the genetic architecture of reproductive traits and disorders across species. In this study, we performed scRNA-seq on Hu sheep ovaries sampled at key developmental stages: embryonic day 90 (E90, onset of primary follicles formation 34 ), young-aged (3 months [M3], pre-puberty; 6 months [M6], post-puberty), and adult (early adulthood; 2 years old [Y2], peak reproductive maturity; Mid-adulthood [Y4], gradual decline in ovarian reserve). Our aim was to investigate the dynamic cellular landscape and transcriptional programs of the developing ovary, providing critical insights into ovarian gene regulation and developmental transitions in sheep. These findings have significant implications for advancing reproductive biotechnology in the sheep industry. Additionally, we conducted cross-species single-cell transcriptomic analyses between humans and sheep, facilitating knowledge transfer between species. By leveraging the conserved transcriptional features of ovarian development, this research reinforces the utility of sheep as a translational model for studying human reproductive traits and diseases.

The authors declare no competing interests.

REAGENT or RESOURCE SOURCE IDENTIFIER Chemicals, peptides, and recombinant proteins 0.25% trypsin Gibco 15050065 0.2% collagenase II Sigma-Aldrich C2-22-1G 40 μm nylon strainer BD Falcon 352340 0.2% FBS Gibco 10091148 Critical commercial assays Chromium Next GEM Single Cell 3′ GEM, Library & Gel Bead Kit v3.1 10xGenomics 1000121 Deposited data Raw ovine scRNAseq data This paper PRJNA1186988 Human ovarian scRNAseq data NCBI GEO: GSE143380 and GSE255690 Human GWAS data Publications Kentistou et al.; Zenin et al.; Liu et al.; Solé-Navais et al.; Day et al.; Rahmioglu et al.; Gallagher et al. 94 , 95 , 96 , 97 , 98 , 99 , 100 Software and algorithms CellRanger https://github.com/10XGenomics/cellranger v7.0.1 Seurat https://satijalab.org/seurat/ v5.0.0 Monocle2 https://cole-trapnell-lab.github.io/monocle-release/ v2.26.0 harmony https://github.com/immunogenomics/harmony v1.1.0 clusterProfiler https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html v4.6.2 ggplot2 https://ggplot2.tidyverse.org/ v3.4.4 pheatmap https://www.rdocumentation.org/packages/pheatmap/versions/1.0.12/topics/pheatmap v1.0.12 pyscenic https://github.com/aertslab/pySCENIC v0.11.2 CellChat https://github.com/sqjin/CellChat v1.6.1 Data analysis pipeline https://github.com/BingruZhao/Singlecell-RNA-seq-of-sheep-ovarian NA Ovary samples were collected from healthy domesticated Hu sheep at critical developmental stages: embryonic day 90 (E90, where primary follicles begin to appear 34 ), young-aged (3 months [M3], pre-puberty; 6 months [M6], post-puberty), and adult stages (2 years old [Y2], peak reproductive maturity; 4 years old [Y4], maintained reproductive function with gradual decline in ovarian reserve). All the sampling experiments were performed in compliance with Basel Declaration agreement. All animal procedures were approved by the Animal Care and Use Committee of the Nanjing Agriculture University (Jiangsu, China) (NO: NJ202306002). All applicable institutional and/or national guidelines for the care and use of animals were followed. All efforts were made to minimize animal suffering. For single-cell isolation, sheep ovaries were finely dissected and enzymatically dissociated using 0.25% trypsin (Gibco, Grand Island, NY, USA) and 0.2% collagenase II (Sigma-Aldrich, C2-22-1G, Darmstadt, Germany) at 37°C for 6 to 8 minutes. The whole ovary at E90 was used to dissociate single cells, and the side of the ovary at other stages was selected with better development, and multiple parts (at least 4 points) of each ovary were mixed into one, and the single cells were dissociated as representative samples at each stage. Following enzymatic dissociation, mechanical dissociation was performed using a pipette to ensure the release of single cells. The resulting cell suspension was filtered through a 40 μm nylon strainer (BD Falcon, 352340, USA) to enrich primary (considered as steady-state) oocytes, centrifuged at 330×g for 10 minutes at 4°C, and resuspended in a base solution containing 0.2% FBS (Gibco). 101 , 102 Post-centrifugation, cells were manually counted three times using the Trypan blue exclusion method to confirm a final concentration of ≥ 2×10 6 cells/ml and cell viability exceeding 85%. The isolated single cells were subsequently processed using a Chromium Controller (10X Genomics) following the manufacturer’s protocol. The synthesis of libraries and RNA sequencing was carried out by Allwegene Technology Co., Ltd (Beijing, China). In brief, cells from each developmental stage were pooled into a single sample and adjusted to a concentration of 1000 cells/μL. Indexed sequencing libraries were prepared using the Chromium Single Cell 3′ Reagent Kits (V2) following the manufacturer’s protocol. The final Single Cell 3′ Libraries included the P5 and P7 primers necessary for Illumina bridge amplification PCR. Barcoded sequencing libraries were quantified using a qPCR assay based on a standard curve (KAPA Biosystems, USA) and further analyzed using an Agilent Bioanalyzer 2100 (Agilent, Loveland, CO, USA). Sequencing of the libraries was performed on an Illumina NovaSeq 6000 platform (Illumina, San Diego, USA) using a custom paired-end sequencing mode of 26 bp (read 1) × 98 bp (read 2). Chromium scRNA-seq raw data were processed for demultiplexing, alignment, and read counting according to the 10XGenomics pipelines using CellRanger v7.0.1 software ( https://github.com/10XGenomics/cellranger ) with default settings. The sequencing FASTQ files were aligned to the sheep reference genome (GCF_016772045.1_ARS-UI_Ramb_v2.0). Gene expression matrices were further analyzed to remove outlier cells and genes with R package Seurat (v5.0.0). 35 We retained cells expressing more than 200 and less than 6,000 genes, a percentage of mitochondrial genes less than 10%, decontX less than 0.9, 103 and log 10 GenesPerUMI greater than 0.8. Genes were retained in the data if they were expressed in ≥3 cells. Furthermore, We used the DoubletFinder package (v2.0.2) 104 to identify and remove potential doublets from the dataset. After applying these quality-control criteria, we used the “merge” function from the Seurat (v5.0.0) 35 package to merge the preprocessed datasets. Additional normalization was performed on the filtered matrix to obtain the normalized count by the “NormalizeData” function. 35 Highly variable genes (n=2000) across single cells were identified using the “FindVariableFeatures” function. 35 Matrix was applied a linear transformation (scaling) using the “ScaleData” function, 35 and principal component analysis was performed to reduce the dimensionality on the top 18 principal components (PCs) which were determined by the “permutationPA” function in the jackstraw R package ( https://github.com/ncchung/jackstraw ). Then, the multiple samples were processed with “harmony” to remove potential batch effects. 105 Cell clustering analysis was then applied by the “FindClusters” function with a resolution of 0.5 and visualized in two dimensions using uniform manifold approximation and projection (UMAP). 35 To accurately define cell types (or subtypes) from ovary tissues, we employed three complementary approaches. First, we automatically annotated the cell clusters using the ScType package. 106 Then, we used the “FindAllMarkers” function from the Seurat (v5.0.0) 35 package to identify positive marker genes with the thresholds of P value 0.25 35 for each cell cluster and analyzed the specific functions of the highly expressed genes for each cell cluster. Thus, we manually annotated these clusters into different cell types (or subtypes) according to the expression of cell type-specific markers as reported in previous literatures (see details in main text), combined with biological functions, and automatically annotation results. Furthermore, differentially expressed genes (DEGs) between each cell group and other groups were identified using the “FindMarkers” function with the thresholds of P value 0.25. Functional enrichment of DEGs for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways was conducted by using ClusterProfiler R package. 107 The similar pipeline was also used in the identification of developmental stage-specifically expressed genes. Single-cell trajectories were analyzed using gene expression matrices across cells by Monocle (v2.20.0). 72 Monocle reduced the space down to one by two dimensions and ordered the cells. Following this ordering, we visualized the cell trajectories within the dimensionally reduced space, which exhibited a tree-like structure with various tips and branches. To explore potential intercellular communication among ovarian cell populations, we used the CellChat package (v1.6.1) 62 to infer ligand-receptor pairs. A ligand or receptor was considered expressed if it showed non-zero expression in at least 20% of cells within a specific cell population. We then constructed a hypothetical cell-cell communication network by linking these expressed ligands to their corresponding receptors across and within major cell populations. The directional signaling from ligand to receptor is depicted through arrows within the chord plot. Additionally, the cumulative number of communication signals transmitted and received by specific cell populations is illustrated by numbered bands in the circular visualization. The circlize R package 108 was used to visualize this hypothetical intercellular communication network. Active TF regulons were predicted using the SCENIC package, 54 employing the standard workflow and the hg38 databases available at https://resources.aertslab.org/cistarget/ . For co-expression analyses, genes were considered if expressed in at least 1% of cells before to the application of the GENIE3 algorithm. 109 To enable cross-species comparative analysis of ovarian single-cell transcriptomic data, we utilized one-to-one orthologous genes between sheep and human genomes. To assess the reproducibility of cell types both within and across species, we employed MetaNeighbor version 1.18. 79 Highly variable genes were identified using the “variableGenes” function, and “MetaNeighborUS” analysis was performed to determine correlations between cell types based on AUROC values. We collected the GWAS summary statistics of 12 human complex traits, including 9 reproductive traits and 3 reproductive diseases 94 , 95 , 96 , 97 , 98 , 99 , 100 with details shown in Table S7 scPagwas 1.3.0, 81 a pathway-based polygenic regression method that integrates scRNA-seq data and GWAS data to identify cell populations critical for complex diseases and traits. To prioritize the top trait-associated genes, we used the scGet_PCC function, which ranks genes based on the Pearson correlation coefficient (PCC) the expression of each gene and the summed genetically associated pathway activity score across all cells. Subsequently, the scPagwas_perform_score function was then applied to analyze pathway activity within each cell type and determine the significance of active pathways using the singular value decomposition (SVD) method. 110 The data analysis and statistical tests employed this this study were primarily done with the R language (version 3.6.5). Student’s t test or two-way analysis of variance (ANOVA) were performed to analyze the data. A p-value of 0.05 was considered significant. Statistical details of the experiments are also provided in method details and Figure legends. No additional resources were generated in the study.

This research was funded by the 10.13039/501100001809 National Natural Science Foundation of China (No. 32441081), the National Key Research and Development Program, China (No. 2021YFF1000702), the Fundamental Research Funds for the Central Universities, China (No. KYT2023002), the 10.13039/501100012453 China Agricultural Research System (CARS-38), the Excellent Postdoctoral Program in Jiangsu Province, China (No. 2023ZB584), and the high-performance computing platform of Bioinformatics Center, 10.13039/501100008562 Nanjing Agricultural University .