A single-cell atlas of the cycling murine ovary

  1. Mary E Morris  Is a corresponding author
  2. Marie-Charlotte Meinsohn
  3. Maeva Chauvin
  4. Hatice D Saatcioglu
  5. Aki Kashiwagi
  6. Natalie A Sicher
  7. Ngoc Nguyen
  8. Selena Yuan
  9. Rhian Stavely
  10. Minsuk Hyun
  11. Patricia K Donahoe
  12. Bernardo L Sabatini
  13. David Pépin  Is a corresponding author
  1. Department of Gynecology and Reproductive Biology, Massachusetts General Hospital, United States
  2. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, United States
  3. Department of Surgery, Harvard Medical School, United States
  4. Howard Hughes Medical Institute, Department of Neurobiology, Harvard Medical School, United States

Abstract

The estrous cycle is regulated by rhythmic endocrine interactions of the nervous and reproductive systems, which coordinate the hormonal and ovulatory functions of the ovary. Folliculogenesis and follicle progression require the orchestrated response of a variety of cell types to allow the maturation of the follicle and its sequela, ovulation, corpus luteum formation, and ovulatory wound repair. Little is known about the cell state dynamics of the ovary during the estrous cycle and the paracrine factors that help coordinate this process. Herein, we used single-cell RNA sequencing to evaluate the transcriptome of >34,000 cells of the adult mouse ovary and describe the transcriptional changes that occur across the normal estrous cycle and other reproductive states to build a comprehensive dynamic atlas of murine ovarian cell types and states.

Editor's evaluation

This manuscript presents an important and useful dataset for understanding cellular and transcriptional dynamics during the estrous cycle in mice. Using single-cell RNA sequencing, the authors' data is compelling, providing new marker genes for different cell types. These data will be useful for understanding ovarian biology and will be of interest to biologists studying other tissues.

https://doi.org/10.7554/eLife.77239.sa0

Introduction

The ovary is composed of a variety of cell types that govern its dynamic functions as both an endocrine organ capable of producing hormones such as sex steroids and a reproductive organ orchestrating the development of follicles, a structure defined by an oocyte surrounded by supporting somatic cells such as granulosa cells and theca cells. Most follicles in the ovary are quiescent primordial follicles, representing the ovarian reserve. Once activated, a primordial follicle grows in size and complexity as it progresses to primary, preantral, and antral stages, adding layers of granulosa and theca cells and forming an antral cavity, until it ultimately ejects the oocyte-cumulus complex at ovulation while the follicular remnants undergo terminal differentiation to form the corpus luteum (CL) (Dunlop and Anderson, 2014). This process necessitates precise coordination of germ cells and several somatic cell types, including granulosa cells, thecal cells, vascular cells, and other stromal cells of the ovary to support the growth of the oocyte until its ovulation or, as is most often the case, undergo follicular atresia. In addition to supporting germ cells, ovarian somatic cells must produce the necessary hormonal cues, as well as coordinate the profound tissue remodeling, necessary to accommodate these dynamic developing structures. For reproductive success to occur, the state of each of these cells must change in a coordinated fashion over the course of the estrous cycle; this allows waves of follicles to grow and mature, ovulation to be triggered precisely, and provides the hormonal support necessary for pregnancy.

Single-cell RNA sequencing (scRNAseq) has been used in a variety of tissues to obtain an in-depth understanding of gene expression and cellular diversity. In the ovary, this technique has allowed us, and others, to explore various physiological processes during early ovarian development and ovarian aging (Zhao et al., 2020; Stévant et al., 2019; Wagner et al., 2020; Niu and Spradling, 2020; Jevitt et al., 2020; Man et al., 2020; Meinsohn et al., 2021; Fan et al., 2019; Wang et al., 2020). For example, Fan et al. cataloged the transcriptomic changes that occur during follicular development and regression and mapped the cell types of the human ovary using surgical specimens (Fan et al., 2019). A primate model has been used to investigate changes in cell types and states that occur in the ovary with aging (Wang et al., 2020). Zhao et al. looked at the formation of the follicle during early embryonic ovarian development to discern the relationship of oocytes to their support cells in formation of follicles (Zhao et al., 2020). We have used scRNAseq to identify inhibitory pathways regulated by anti-Müllerian hormone (AMH) during the first wave of follicular growth in the murine ovary (Meinsohn et al., 2021). While all these studies have helped establish a static framework to understand the major cell types in the ovary, they fail to describe the dynamic nature of cell states across the reproductive cycle, known as the estrous cycle. The estrous cycle in mice is analogous to the human menstrual cycle, which both reflect follicle development in the ovary. In mice, this cycle lasts 4–5 days and is composed of four different phases known as proestrus, estrus, metestrus, and diestrus. The murine proestrus is analogous to the human follicular stage and leads to ovulation at estrus. Metestrus and diestrus are analogous to early and late secretory stages of the reproductive cycle in humans, which are orchestrated by production of progesterone by the CL (Ajayi and Akhigbe, 2020).

To understand more fully the dynamic effects of cyclic endocrine, autocrine, and paracrine signals on ovarian cell states, we performed high-throughput scRNAseq of ovaries from adult mice across a physiological spectrum of reproductive states. Ovaries were harvested from mice in the four phases of the normal estrous cycle: proestrus, estrus, metestrus, and diestrus. Additionally, ovaries were evaluated from mice that were either lactating or non-lactating 10 days post-partum, and from randomly cycling adult mice to increase the diversity of cell states represented in the dataset. Herein, we (1) describe the previously unrecognized complexity in the ovarian cellular subtypes and their cyclic expression states during the estrous cycle, and (2) identify secreted factors that cycle and thus could represent potential biomarkers for staging.

Results

scRNA-seq of adult mouse ovaries across reproductive states

To survey the dynamic transcriptional landscape of ovaries at the single-cell level across a range of physiological reproductive states in sexually mature female mice, we isolated the ovaries (four mice per group) at each stage of estrous cycling (proestrus, estrus, metestrus, and diestrus), post-partum non-lactating (PPNL) (day 10 post-partum, with pups removed on the day they were born), post-partum lactating (day 10 post-partum, actively lactating with pups), and non-monitored adult mice to increase sample diversity and cell counts. Following enzymatic digestion of the ovaries, we generated single-cell suspensions and sorted them by microfluidics using the inDROP methodology (Klein et al., 2015), targeting 1500 cells per animal. Resulting libraries were indexed and combined for sequencing (Figure 1A).

Figure 1 with 1 supplement see all
Single-cell RNA sequencing of cycling mouse ovaries.

(A) Schematic of the single-cell sequencing pipeline. (B) Uniform manifold approximation and projection (UMAP) plot featuring the different clusters of the ovary and their composition by stage of the estrous cycle, lactating status, or unmonitored. (C) Heatmap of the top 10 markers of each cluster by fold change.

Following dimensionality reduction and clustering using the Seurat algorithm, we identified multiple clusters which could be combined to represent the major cell categories of the ovary (Figure 1B). To assign cell type identity, we used cluster-specific markers which were previously described in other studies or newly identified makers later validated by RNA in situ (Supplementary file 2). The largest groups of clusters consisted of granulosa cells (N=17627 cells) and mesenchymal cells of the ovarian stroma (N=10825 cells). Other minor cell types were identified including endothelial cells (N=3501 cells), ovarian surface epithelial cells (N=1088 cells), immune cells (N=1649 cells), and oocytes (N=22 cells), altogether recapitulating all the major cell types of the ovary (Figure 1—figure supplement 1A). Oocytes were poorly represented in the dataset due to cell size limitations of inDROP, likely restricting our sampling to small oocytes of primordial follicles (Figure 1B). To characterize more fully the transcriptional signatures of the identified cell types, we evaluated a heatmap of marker gene expression across the major categories of cell types and states (Figure 1C). Cells were also classified depending on the stage of the estrous cycle or lactating states in which the ovaries were collected (Figure 1B). Morphological differences between the stages of proestrus, estrus, metestrus, diestrus, and also post-partum lactating and non-lactating, were documented in Figure 1—figure supplement 1B. The granulosa, mesenchyme, and epithelium clusters were isolated and reanalyzed to identify subclusters.

Single-cell sequencing reveals heterogeneity within granulosa and mesenchymal cell clusters

Cellular diversity of mesenchymal cells

The mesenchymal cluster was the second largest cluster identified in our analysis. Based on prior studies and conserved marker expression (Fan et al., 2019; Wang et al., 2020), we were able to identify subclusters within mesenchymal cells and their relative abundance (percentage) as follows: early theca (16.8%), which formed the theca interna of preantral follicles; steroidogenic theca (13.2%), which formed the theca interna of antral follicles; smooth muscle cells (10.2%), which were part of the theca externa of both antral and preantral follicles; pericytes (6.2%), which surrounded the vasculature; and two interstitial stromal cell clusters, one composed of steroidogenic cells (28.7%) and the other of fibroblast-like cells (24.9%), which together constituted the bulk of the ovarian volume (outside of follicles). These subclusters can be seen in Figure 2A, with the top five expressed markers of each subcluster described in the Figure 2B heatmap and the top 10 listed in Supplementary file 3.

Figure 2 with 1 supplement see all
Identification of the different cell types of the mesenchyme cluster.

(A) UMAP plot featuring the different cell subclusters belonging to the mesenchyme cluster. (B) Heatmap of the top five markers of each subcluster by fold change. (C) Validation of the identity of mesenchyme subcluster by UMAP-plots (cluster of interest circled) and RNA in situ hybridization.

Distinct transcriptional signatures were identified in each of these mesenchymal subclusters (Figure 2B); to confirm the presumed identity and histology of these cell types (detailed in Figure 1—figure supplement 1A), we validated markers prioritized by highest fold-change expression, highest differential percent expression, and lowest p value (Figure 2C).

For the theca interna, the two clusters identified reflected the stage of development of the follicle: early thecal cells could be defined by their expression of hedgehog-interacting protein (Hhip) and were histologically associated with preantral follicles. Meanwhile, the steroidogenic theca cells were identified by their expression of cytochrome P450 family 17 subfamily A member 1 (Cyp17a1), an essential enzyme for androgen biosynthesis (Richards et al., 2018); they were found in antral follicles (Figure 2C). The theca externa is a connective tissue rich in extracellular matrix situated on the outermost layer of the follicle (Figure 1—figure supplement 1A), containing fibroblasts, macrophages, blood vessels, and abundant smooth muscle cells, which we identified based on expression of microfibril-associated protein 5 (Mfap5) by RNA in situ hybridization (Figure 2C). To validate the identity and histology of these smooth muscle cells, we performed RNAish/IHC colocalization of Mfap5 and actin alpha 2 (Acta2), another marker of smooth muscle, which confirmed their position within the theca externa. In contrast, Hhip, which was expressed in theca interna (both immature and steroidogenic), did not colocalize with Acta2 (Figure 2—figure supplement 1A-C). These results suggest Mfap5 labels smooth muscle cells of the theca externa more specifically than Acta2; these cells are thought to perform a contractile function during ovulation (Young and McNeilly, 2010).

Lastly, the bulk of the ovarian interstitial stromal space was made up of two closely related cell types which could not be differentiated by specific dichotomous markers but rather were distinguished based on relative expression of ectonucleotide pyrophosphatase/phosphoiestrase 2 (Enpp2) (Figure 2C). While Enpp2+ cells represented fibroblast-like stromal cell, Enpp2− interstitial cells were enriched for expression of genes such as Patch1 (Ptch1), a member of the hedgehog-signaling pathway, an important regulator of ovarian steroidogenesis (Spicer et al., 2009), suggesting these represented steroidogenic stromal cells. Indeed, the steroidogenic activity of this stromal cell cluster was further confirmed by its high relative expression of other genes associated with steroidogenesis including cytochrome P450 family 11 subfamily A member 1 (Cyp11a1), hydroxy-delta-5-steroid dehydrogenase, 3 beta- and steroid delta-isomerase 1 (Hsd3b1), cytochrome P450 family 17 subfamily A member 1 (Cyp17a1), steroid 5 alpha-reductase 1 (Srd5a1), along with other markers such as potassium two pore domain channel subfamily K member 2 (Kcnk2) (Figure 2—figure supplement 1E, F). In contrast the fibroblast-like stromal cluster had enriched expression of many extracellular matrix genes such as collagen type I alpha 1 chain (Col1a1), collagen type V alpha 1 chain (Col5a1), Lumican (Lum), lysyl oxidase like 1(Loxl1) (Muhl et al., 2020) along with other known markers of fibroblasts such as C-X-C motif chemokine ligand 14 (Cxcl14) (Lu et al., 2016) and WT1 transcription factor (Wt1) (He et al., 2008; Figure 2—figure supplement 1E, F).

The identities of above-described mesenchymal clusters matched those of known ovarian stromal cell types based on expression of previously reported markers such as Desmin (Des) for pericytes (Hughes and Chan-Ling, 2004), steroidogenic acute regulatory protein (Star) for the steroidogenic theca (Kiriakidou et al., 1996), cellular communication network factor 1 (Ccn1) for smooth muscle cells (Yang et al., 2018b), receptor activity-modifying protein 2 (Ramp2) for the early theca cells surrounding preantral follicles (Hatzirodos et al., 2015), and finally C-X-C motif chemokine ligand 12 (Cxcl12) which we found in both stromal clusters (Porcile et al., 2005) as illustrated in Figure 2—figure supplement 1D.

Cellular diversity of granulosa cells

To explore further the cellular heterogeneity within developing follicles (listed in Figure 1—figure supplement 1A), we investigated the subclustering of granulosa cells based on their transcriptional profile. Consistent with previous reports, we could distinguish discrete granulosa cell states in follicles based on their stage of development (Zhao et al., 2020; Fan et al., 2019; Gallardo et al., 2007). Granulosa cells could be subdivided into eight main categories: preantral-cumulus (27.3%), antral-mural (21.8%), luteinizing mural (4.8%), atretic (22.6%), mitotic (14.4%), regressing CL (3.7%), and active CL (5.4%) (Figure 3A). Supplementary file 4 lists the top 10 markers for each of these clusters. Distinctive gene expression programs were identified in the granulosa cell subclusters, as visualized in the heatmap (Figure 3B), from which we selected potential markers for validation.

Figure 3 with 1 supplement see all
Identification of the different cell types in the granulosa cluster.

(A) UMAP plot featuring the different cell subclusters belonging to the granulosa cluster (specific subcluster circled in each UMAP). (B) Heatmap of the top five markers of each cluster by fold change. (C) Validation of the identity of granulosa subclusters by UMAP-plots and RNA in situ hybridization.

Early preantral granulosa cells, and those constituting the cumulus oophorus of antral follicles, could be identified by their shared expression of markers such as potassium channel tetramerization domain (Kctd14) (Figure 3C), which we had previously shown to be expressed by preantral follicles (Meinsohn et al., 2021). In contrast, mural granulosa cells of antral follicles expressed distinct markers (Supplementary file 4) such as male-specific transcription in the developing reproductive organs (Mro) (Figure 3C). Luteinizing mural granulosa cells could be identified by the expression of previously established markers (Supplementary file 4) and oxytocin receptor gene (Oxtr) which we propose as a highly specific marker for this cell type, a likely target of the surge in oxytocin during estrus (Ho and Lee, 1992; Figure 3C). Furthermore, we identified two different clusters that we hypothesize represent cell states of the CL, either active or regressing, which both expressed nuclear paraspeckle assembly transcript 1 (Neat1), a known marker of CLs (Nakagawa et al., 2014). To confirm the active and regressing CL cell states, we investigated the expression of Top2a, a mitotic marker (Donadeu et al., 2014), which was enriched in the active CL cluster, and Cdkn1a, a cell cycle exit and senescence marker (Ock et al., 2020), which was enriched in the regressing cluster (Figure 3—figure supplement 1B, C). Moreover, when examining the composition of clusters depending on the reproductive stage, the regressing CL cluster was found to be composed mostly of cells derived from the Postpartum non lactating (PPNL) samples (Figure 3—figure supplement 1E), which overexpressed markers related to CL regression (Talbott et al., 2017; Figure 3—figure supplement 1F), consistent with a post-partum effect of prolactin. Finally, two relatively abundant granulosa cell states could be identified based on marker expression: mitotic granulosa cells could be found in both preantral and antral follicles and were defined by their expression of Top2a, and atretic granulosa cells, which expressed markers consistent with follicular atresia and apoptosis such as phosphoinositide-3-kinase-interacting protein 1 (Pik3ip1), nuclear protein 1, transcriptional regulator (Nupr1), growth arrest and DNA damage inducible alpha (Gadd45a), vesicle amine transport 1 (Vat1), transgelin (Tagln), and melanocyte-inducing transcription factor (Mitf) (Terenina et al., 2017; Figure 3C, Figure 3—figure supplement 1A, Supplementary file 4). Furthermore, we propose growth hormone receptor (Ghr), which was highly specific to this cluster, as a specific marker of atretic follicles, which warrants further investigation of the role of growth hormone in this process (Figure 3C).

Cellular states in the ovarian surface epithelium

The epithelial cluster was composed of 1088 ovarian surface epithelium (OSE) cells, which could be further subdivided into two clusters (Figure 4A): the larger one composed of non-dividing epithelium cells (96%), and a smaller cluster (4%), composed of mitotic epithelium. The latter was characterized by proliferation markers such as thymidine kinase 1 (Tk1) (Liu et al., 2019), Rac GTPase-activating protein 1 (Racgap1) (Yang et al., 2018a), Top2a, casein kinase 1 (Ck1) (Gao et al., 2021), protein regulator of cytokinesis 1 (Prc1) (Liang et al., 2019), ubiquitin-conjugating enzyme E2 C (Ube2c) (Xiong et al., 2019), and baculoviral IAP repeat containing 5 (Birc5) (Xu et al., 2021; Figure 4A and B). Interestingly, the proliferating subcluster of OSE was almost exclusively composed of cells from the estrous stage (Figure 4C and D), consistent with their transient amplification during ovulatory wound closure (Mara et al., 2020).

Identification of epithelial subclusters.

(A) UMAP plot of the surface epithelium cluster showing two subclusters: epithelium and mitotic epithelium (circled in black). (B) Heatmap of proliferation markers expressed in the proliferating epithelium cluster. (C) UMAP plot of the cellular composition of the epithelium subclusters by reproductive state (mitotic subcluster circled in the estrous state). (D) Expression of proliferation markers depending on the phase of the estrous cycle.

Granulosa cell transcriptome is most dynamic during the proestrus/estrus transition

To identify changes in cell states associated with the stages of the estrous cycle, we focused on the granulosa cell subclusters, given the importance of follicular maturation in coordinating this process (illustrated in Figure 1—figure supplement 1B). When comparing the composition of granulosa cell subclusters by estrous stage, we found that some clusters were dominated by cells from either the proestrous or estrous samples, particularly the clusters corresponding to ‘antral/mural’ and ‘periovulatory’ clusters, respectively (Figure 5A and B). A volcano plot analysis confirmed that the transition between these two stages was characterized by 24 significantly upregulated and 10 significantly downregulated markers (Figure 5C), which together with the transition from estrus to metestrus represents the largest change in gene expression. In contrast, few genes were found to change significantly during the transition from metestrus to diestrus, or diestrus to proestrus (Figure 5—figure supplement 1A). Gene ontology analysis revealed that the most significantly differentially regulated pathways between the proestrous and estrous phases were related to ovarian matrix remodeling and steroidogenesis and hormones production (Figure 5—figure supplement 1B). To validate the genes with significant changes in expression identified within the single-cell sequencing dataset, we performed quantitative PCR (qPCR) on whole-ovary samples at the proestrus to estrus transition, including the steroid biosynthesis markers cytochrome P450 family 19 subfamily A member 1 (Cyp19a1, p=0.0029, proestrus to estrus), Star protein (p=0.0187, proestrus to estrus), serum- and glucocorticoid-inducible kinase-1 (Sgk1, p=0.0056, proestrus to metestrus), as well as matrix remodeling genes such as regulator of cell cycle (Rgcc, p=0.0441, proestrus to estrus), tribbles pseudokinase 2 (Trib2, p=0.0023, proestrus to estrus) (Figure 5D and E), and immediate early genes, fos proto-oncogene (Fos), jun proto-oncogene (Jun, p=0.0022, proestrus to estrus), jun proto-oncogene B (Junb, p=0.0069, proestrus to diestrus), and early growth response 1 (Egr1, p=0.0504 estrus to diestrus), which represent a family of genes thought to be involved in wound repair, a sequela of ovulation (Florin et al., 2006; Wu et al., 2009; Martin and Nobes, 1992; Yue et al., 2020; Figure 5—figure supplement 1C). Transcriptional gene expression changes were found to be concordant between the scRNAseq data and whole-ovary transcripts quantified by qPCR.

Figure 5 with 1 supplement see all
Gene expression in granulosa cells by estrous stage.

(A) UMAP plot featuring estrous cycle stages in the granulosa cell cluster. (B) UMAP plots featuring each of the estrous cycle phases individually (proestrus and estrus enriched subclusters circled in black). (C) Volcano plot of genes differentially expressed between proestrus and estrous stages. (D) DotPlot of differentially expressed markers between proestrus and estrus. (E) Quantitative PCR (qPCR) validation of differentially expressed genes involved in extracellular matrix remodeling and steroidogenesis markers (n=5 per group, mean ± SEM, *p<0.05, **p<0.01, ***p<0.005, and ****p<0.001).

Identification and validation of secreted biomarkers varying throughout the estrous cycle

To identify new biomarkers that vary as a function of the estrous cycle and that could be used for staging in reproductive medicine, we screened for differentially expressed secreted factors (DAVID Bioinformatics Resources) (Sherman et al., 2022; Huang et al., 2009), which would therefore be potentially measurable in the blood. Furthermore, to ensure specificity, we prioritized genes expressed specifically in the granulosa or ovarian mesenchymal clusters and not highly expressed in other tissues based on their GTEX profile (GTEx Consortium, 2013; Supplementary file 5). As a primary screen, we first validated our ability to detect gene expression changes by estrous stage using whole-ovary qPCR analysis in a separate set of staged mice (N=4 per group). Whole-ovary qPCR successfully detected expression changes of estrous cycle markers such as luteinizing hormone/choriogonadotropin receptor (Toms et al., 2017) (Lhcgr, p=0.0281 estrus to metestrus) and progesterone receptor (Pgr, p=0.0096, proestrus to estrus) (Kubota et al., 2016; Figure 6B). Using this method, we validated a set of significantly upregulated secreted markers in the proestrous to estrous transition, most prominent of which were natriuretic peptide C (Nppc, p=0.0022 proestrus to estrus) and inhibin subunit beta-A (Inhba, p=0.0067, proestrus to estrus) (Figure 6A and B). Similarly, tubulointerstitial nephritis antigen like 1 (Tinagl1) and serine protease 35 (Prss35) were secreted markers significantly upregulated in estrus compared to their level of transcription in proestrus in the scRNAseq dataset (Figure 6A) and by qPCR (Tinagl1, p=0.0081, proestrus to estrus; Prss35, p=0.0008, proestrus to estrus) (Figure 6B). In situ RNA hybridization showed that, as expected, these markers were mostly expressed in mural granulosa cells of antral follicles, while Nppc was expressed in both mural and cumulus cells (Figure 6C).

Identification and validation of new secreted estrous staging markers.

(A) Expression of granulosa cell transcripts varying by estrous cycle stage. (B) Validation of significantly up- and downregulated transcripts of secreted estrous staging markers by qPCR (n=5 per group, mean ± SEM, *p<0.05, **p<0.01, and ***p<0.005). (C) Localization of estrous staging markers by in situ hybridization (RNAscope) in ovarian sections. (D) Quantification of circulating estrous staging markers proteins in the blood by enzyme-Linked Immunosorbent assay (ELISA) (n=5 per group, mean ± SEM, *p<0.05, and ***p<0.005). (E) Summary of the timing of expression of estrous staging markers in the blood.

To evaluate the feasibility of measuring the secreted PRSS35, NPPC, TINAGL1, and activin A proteins in the serum for staging, we performed ELISAs in mice at each stage of the estrous cycle (Figure 6D). We found that the activin A concentration in the serum was significantly increased between the diestrous and proestrous stages (p=0.0312) and peaked at the proestrous stage (Figure 6D). The Inhba transcript, which encodes for the activin and inhibin beta-A subunit, had a similar temporal expression profile (Figure 6B). Circulating PRSS35 levels were lowest at the metestrous stage and were significantly increased during the transition to diestrus (p=0.0009) and remained significantly elevated until the proestrus (Figure 6D). In contrast, the Prss35 transcript was significantly induced earlier at estrus (Figure 6B). The serum concentrations of TINAGL1, which was lowest at the diestrous and proestrous stages, was significantly increased during the transition between proestrus and metestrus, peaking in estrus (p=0.0142) (Figure 6D). This temporal pattern of expression was recapitulated at the transcriptional level by qPCR and scRNAseq (Figure 6A and B). Finally, we observed a trend for serum protein concentrations of NPPC to be lowest at the proestrous and estrous stages and increase during the metestrous and diestrous stages (Figure 6D), although the differences were not statistically significant (p=0.0889, estrus to metestrus). Importantly, these data provide a proof of concept that four markers could be used to monitor estrous cycle progression when measured in conjunction in the blood (Figure 6E).

Discussion

scRNAseq has been used to catalog the transcriptomes of a variety of tissues in several species, across different physiological states (Hwang et al., 2018). Herein, we used scRNAseq to survey the cellular diversity and the dynamic cell states of the mouse ovary across the estrous cycle and other reproductive states such as post-partum lactating (PPL) and post-partum non-lactating (PPNL).

The most significant changes in composition and cell states were identified in granulosa cells, particularly as they cycled through the estrous stages, reflective of their important role in cyclic follicular maturation and hormone production. Early preantral follicle numbers are thought to be relatively stable across the estrous cycle (Deb et al., 2013), given that they are largely unresponsive to gonadotropins (Richards, 1980), in contrast to antral follicles, whose numbers and size are more variable (Deb et al., 2013). Indeed, while subclusters such as ‘preantral granulosa cells’ were equally represented in samples from proestrus, metestrus, and diestrus, others, such as the ‘luteinizing mural’ cluster, were dominated by cells derived from one stage (in this case ‘estrus’). Genes enriched in this cluster had been previously reported to be involved in the ovulatory process and regulated by the luteinizing hormone (LH) surge, including markers of terminal differentiation and steroidogenesis such as Smarca1 (Lazzaro et al., 2006), Cyp11a1 (Irving-Rodgers et al., 2009), metallothionein 1 (Mt1), and metallothionein 2 (Mt2) (Wang et al., 2018; Supplementary file 4). Other genes enriched in this subcluster include Prss35 (Wahlberg et al., 2008) and Adamts1 (Lussier et al., 2017; Sayasith et al., 2013), which had previously been identified as playing a role in the follicular rupture necessary for ovulation.

Interestingly, we found that granulosa cells of preantral follicles and cumulus cells of antral follicles clustered together and shared markers that distinguished them from mural granulosa cells. For example, Kctd14, a member of the potassium channel tetramerisation domain-containing family, was expressed in granulosa cells during the initial growth of early preantral follicles, but also specifically expressed only in cumulus cells, but not mural granulosa cells, of larger antral follicles. Intriguingly, we have previously shown that AMH (anti- Müllerian hormone, a.k.a Müllerian Inhibiting Substance), which is specifically expressed by cumulus cells (Diaz et al., 2007) in antral follicles, regulates KCTD14 expression in preantral follicles (Meinsohn et al., 2021). This conservation of cellular state and marker expression from preantral granulosa cells to cumulus cells of antral follicles suggests a continuous lineage, potentially defined and maintained by the close interaction with the oocyte (Diaz et al., 2007). This interpretation is consistent with the presence of a differentiation fork in the granulosa cell lineage during antrum formation, which would give rise to a distinct mural granulosa cell fate poised to respond to the LH surge. Indeed, the periovulatory granulosa cell state was identified based on its expression of genes regulated by LH such as Smarca1 (Lazzaro et al., 2006), and we propose Oxtr as a specific marker for these cells (Figure 3—figure supplement 1D). Oxtr expression was found only in the mural cells of large Graafian follicles, suggesting it indeed corresponds to an LH-stimulated mural granulosa cell state.

After ovulation, these LH-stimulated mural granulosa cells, along with the steroidogenic theca cells, terminally differentiate into luteal cells and form the corpus luteum (CL). The CL is a transient structure with highly active steroid biosynthesis, providing the progesterone to maintain pregnancy (Duncan, 2021). In absence of implantation, the CL degenerates (Noguchi et al., 2017). We found this progression of the CL to be recapitulated at the transcriptional level, leading to two luteal subclusters: active CL and regressing CL. While the active CL was characterized by expression of proliferation markers (Top2a) in addition to steroidogenic enzymes, the regressing CL expressed the cell cycle inhibitor and senescence maker Cdkn1a, along with luteolysis markers such as syndecan 4 (Sdc4), claudin domain containing 1 (Cldnd1), and BTG anti-proliferation factor 1 (Btg1) (Talbott et al., 2017; Zhu et al., 2013). The distinct expression signatures observed in these two clusters may provide insights into the molecular basis of luteolysis and warrant further investigation.

The mesenchymal cluster was also surprisingly complex and variable across the estrous cycle, reflecting stromal remodeling and other physiological functions supporting follicle growth, steroid hormone production, and ovulation. For example, during follicle maturation, the ovarian stromal cells adjacent to the developing follicle differentiates into the theca, which is ultimately responsible for steroid hormone biosynthesis and therefore underlies the cyclic hormone production of the ovary (Ryan and Petro, 1966). Herein, we identified two thecal clusters, designated as early and steroidogenic theca. Early theca was defined by markers such as Hhip (Richards et al., 2018; Hummitzsch et al., 2019), mesoderm-specific transcript (Mest) (Fan et al., 2019), and patched 1 (Ptch1) (Fan et al., 2019; Richards et al., 2018) and given its association with small follicles, presumed to be immature and the precursor to steroidogenic theca. As the follicle matures and the antrum forms, this layer becomes more vascularized and differentiates into theca interna, which is steroidogenic. This steroidogenic theca cluster was readily identifiable through its expression of steroidogenic enzymes such as Hsd3b1, Cyp17a1, Cyp11a1, and also well-established markers such as ferredoxin-1 (Fdx1) and prolactin receptor (Prlr) (Fan et al., 2019; Grosdemouge et al., 2003). Interestingly, we confirmed the presence of steroidogenic interstitial stromal cells which also likely contribute to sex steroid production in the ovary. Indeed, such cells likely represent the precursors of the theca interna (Sheng et al., 2022; Kinnear et al., 2020). Smooth muscle cells, which are part of the theca externa, were identified by their expression of structural proteins such as Mfap5, myosin heavy chain 11 (Myh11), Tagln, and smooth muscle actin (Sma or Acta2) (Zhao et al., 2020). In contrast to mice, human smooth muscle cells are thought to express high levels of collagen (Fan et al., 2019), which we did not observe here. Another species difference between mice and humans was the expression of aldehyde dehydrogenase 1 family member A1 (Aldh1a1), which we found primarily in the steroidogenic theca cluster, while it is presumably enriched in the theca externa in humans (Fan et al., 2019).

Ovulation is associated with a dramatic remodeling of the ovary, including the subsequent ovulatory wound repair. We identified fibroblast-like cells in the ovarian stroma expressing many of the extracellular matrix protein known to play a role in these processes (Mara et al., 2020; Duffy et al., 2019). Another important player in ovulatory would repair is the ovarian surface epithelium (OSE), a simple mesothelial cell layer that covers the surface of the ovary and must dynamically expand to cover wound (Hartanti et al., 2019; Xu et al., 2011). The OSE cluster could be identified based on well-established markers such as keratin (Krt) 7, 8, and 18 (Kenngott et al., 2014) and was represented by only 3% of all cells in our dataset, which could be further subdivided in proliferative and non-proliferative states. As expected from their function in ovulatory wound repair, dividing OSE was enriched during estrus. Furthermore, genes associated with wound healing such as galectin 1 (Lgals1) (Lin et al., 2015) were also significantly upregulated in estrus. Similarly, the expression of the immediate-early genes Fos, Jun, Junb, and Egr1 was variable during the estrous cycle, following a common pattern of strong downregulation at estrous compared to the other stages, consistent with their temporal expression during the repair of other tissues such as the cornea (Okada et al., 1996).

Finally, to take advantage of this rich dataset, we sought to identify secreted markers which varied in abundance during the estrous cycle and could thus be used as staging biomarkers in assisted reproduction. We identified and prioritized four secreted biomarkers, expressed in mouse but also human ovaries, which varied significantly during different transitions of the estrous cycle, namely Inhba (Wijayarathna and de Kretser, 2016), Prss35 (Wahlberg et al., 2008; Li et al., 2015), Nppc (Zhang et al., 2010; Xi et al., 2021), and Tinagl1 (Akaiwa et al., 2020; Kim et al., 2010).

Activin A is a secreted protein homodimer translated from the Inhba transcript that is a crucial modulator of diverse ovarian functions including pituitary feedback, and whose expression level depends highly on the stage of the estrous cycle (Chang and Leung, 2018). Quantification of activin A protein in the serum by ELISA revealed elevated levels in the blood during both proestrus and estrus, which is consistent with studies of other species such as ewes (OConnell et al., 2016). Importantly, the protein product of Inhba, the inhibin beta-A subunit, can be incorporated into other protein dimers, such as activin BA, and inhibin A, which were not measured in this study and may also represent cycling biomarkers.

The serine protease 35 transcript was expressed in the theca layers of preantral follicles and induced in granulosa cells of preovulatory follicles and all stages of the corpora lutea, peaking at the estrous stage according to qPCR, leading us to speculate that it may be involved tissue remodeling during ovulation and CL formation (Wahlberg et al., 2008). In contrast, the PRSS35 protein levels were highest in the diestrus and proestrus stages as determined by ELISA, suggesting other tissue sources of PRSS35 or an offset in peak protein levels due to delays in accumulation of the protein in the circulation.

The natriuretic peptide precursor C (NPPC) protein is a peptide hormone encoded by the Nppc gene. Nppc has been reported to be expressed by mural granulosa cells, while its receptor Npr2 is expressed by cumulus cells (Zhang et al., 2010). The pair acts on developing follicles by increasing the production of intracellular cyclic guanosine monophosphate and maintains oocyte meiotic arrest during maturation. Upon downregulation of this pathway, the oocyte can escape meiotic arrest and ovulate (Celik et al., 2019). This close relationship with the ovulatory process makes Nppc an attractive marker to predict ovulation. Herein, qPCR analysis revealed that Nppc was highest in the ovary at proestrus and was quickly and significantly downregulated at estrous, probably in response to the increased levels of LH which in turn inhibit the Nppc/Npr2 system (Celik et al., 2015). In contrast, there was a trend for the circulating NPCC peptide to be highest in metestrus and diestrus, albeit not in a statistically significant way.

Finally, we evaluated the level of transcription and protein expression of the matricellular factor Tinagl1. We found both the Tinagl1 transcript and the circulating TINAGL1 protein in the blood to be highest during estrous, thus coinciding with ovulation, with a pattern of expression consistent with expression by mural granulosa cells of antral follicles. While the role of TINAGL1 in the ovary has not been extensively investigated, it has been associated with delayed ovarian collagen deposition and increased ovulation in aging Tinagl1 knock-out mice (Akaiwa et al., 2020).

Those four potential cyclic biomarkers, activin A, PRSS35, NPPC, and TINAGL1, provide a proof of concept that a deeper understanding of transcriptional changes at the single-cell level may translate into useful applications in assisted reproduction. It will be of interest to follow up the findings of cyclic expression of these four markers, particularly in combination as an index, for the purpose of staging and predicting ovulation timing in humans and other species .

In summary, this study outlines the dynamic transcriptome of murine ovaries at the single-cell level and across the estrous cycle and other reproductive states, and extends our understanding of the diversity of cell types in the adult ovary. We identified herein novel biomarkers of the estrous cycle that can be readily measured in the blood and may have utility in predicting staging for assisted reproduction. This rich dataset and extensive validation of new molecular markers of cell types of the ovary will provide a hypothesis-generating framework of dynamic cell states across the cycle with which to elucidate the complex cellular interactions that are required for ovarian homeostasis.

Materials and methods

Key resources table
Reagent type (species) or resourceDesignationSource or referenceIdentifiersAdditional information
Genetic reagent (Mus musculus)C57BL/6-Tg(UBC-GFP)30Scha/JJackson Laboratorystock #004353
AntibodySmooth muscle alpha action (SMA) (Rabbit polyclonal)Abcam#5694Dilution: 1:300
Commercial assay or kitACTIVIN A commercial ELISARnD systems#DAC00B
Commercial assay or kitNPPC commercial ELISANovus Bio#NBP2-75790
Commercial assay or kitTinagl1 commercial ELISALS-Bio#LS-F49684
Commercial assay or kitPRSS35 commercial ELISAMybiosource#MBS9717242
Commercial assay or kitRNA scope 2.5 HD Duplex detection kitACD bio#322500
Commercial assay or kitRNA scope 2.5 HD red detection kitACD bio#322360
Commercial assay or kitThe target retrieval and protease plus reagentsACD bio#322330
OtherCdkn1a (M. musculus) NM_007669.4ACD bio# 408551RNAscope probe
OtherCxcl14 (M. musculus) NM_019568.2ACD bio#459741RNAscope probe
OtherCyp11a1 (M. musculus) NM_019779.4ACD bio#809181RNAscope probe
OtherCyp17a1 (M. musculus) NM_007809.3ACD bio#522611RNAscope probe
OtherGhr (M. musculus) NM_010284.3ACD bio#464951RNAscope probe
OtherHhip (M. musculus) NM_020259.4ACD bio#448441RNAscope probe
OtherInhba (M. musculus) NM_008380.1ACD bio# 455871RNAscope probe
OtherKcnk2 (M. musculus) NM_001159850.1ACD bio#440421RNAscope probe
OtherKctd14 (M. musculus) NM_001136235.1ACD bio#517811RNAscope probe
OtherMfap5 (M. musculus) NM_015776.2ACD bio#490211RNAscope probe
OtherMro (M. musculus) NM_001305882.1ACD bio#491181RNAscope probe
OtherNeat1 (M. musculus) NR_003513.2ACD bio#440351RNAscope probe
OtherNppc (M. musculus) NM_010933.5ACD bio# 493291RNAscope probe
OtherOnecut2 (M. musculus) NM_194268.2ACD bio#520541RNAscope probe
OtherOxtr (M. musculus) NM_001081147.1ACD bio#412171RNAscope probe
OtherPrss35 (M. musculus) NM_178738.3ACD bio#492611RNAscope probe
OtherRunx1 (M. musculus) NM_001111021.1ACD bio#406671RNAscope probe
OtherTinagl1 (M. musculus) NM_001168333ACD bio#312621RNAscope probe
OtherTop2a (M. musculus) NM_011623.2ACD bio# 491221RNAscope probe
OtherWt1 (M. musculus) NM_144783.2ACD bio#432711RNAscope probe
Software and algorithmR version 4.1.3R Project for Statistical Computinghttps://scicrunch.org/resolver/SCR_001905
Software and algorithmSeurat package 4.1.0R toolkit for single-cell genomicshttps://satijalab.org/seurat/articles/install.html
Software and algorithmBZ-X800 analysis softwareKeyencehttps://www.keyence.com/landing/microscope/lp_fluorescence.jsp
Software and algorithmGraphPad Prism, version 9.2.0Graphpad

Mice

Animal experiments were carried out in 6–8 weeks old C57BL/6 mice obtained from Charles River Laboratory, approved by the National Institute of Health and Harvard Medical School Institutional Animal Care and Use Committee, and performed in accordance with experimental protocols 2009N000033 and 2014N000275 approved by the Massachusetts General Hospital Institutional Animal Care and Use Committee.

For the analysis of transcriptional changes in ovaries of cycling mice, animals were housed in standard conditions (12/12 hr light/dark non-inverting cycle with food and water ad libitum) in groups of five females with added bedding from a cage that previously housed an adult male mouse to encourage cycling. Estrous stage was determined by observation of the vaginal opening and by vaginal swabs done at the same time daily, as previously described (Kano et al., 2017). Each mouse was monitored for a minimum of 2 weeks to ensure its cyclicity. Four mice were sacrificed in each of the four phases of the estrous cycle and labeled as being from experimental batch ‘cycling’. An additional eight mice were included in the analysis and labeled as being from experimental batch ‘lactating’. Four of these mice were lactating at day 10 post-partum, and four were 10 days post-partum with pups removed at delivery. Four additional mice were not monitored for cycling and included to increase sample size and diversity.

Additional mice were monitored throughout the estrous cycle to collect ovaries at each stage (groups of N=5 for proestrus, estrus, metestrus, and diestrus) for gene validation. Paired ovaries were collected from each staged mouse: one was used to extract mRNA for qPCR, while the other was fixed in 4% paraformaldehyde for RNAish (RNAscope) or immunohistochemistry to validate gene expression.

Superovulation

Request a detailed protocol

To stimulate superovulation, mature female mice (6–9 weeks C57BL/6) were injected intraperitoneally (IP) with 5 IU of pregnant mare serum gonadotropin (PMSG; Calbiochem, San Diego, CA, USA), followed 48 hr later by 5 IU of human chorionic gonadotropin (hCG; Millipore Sigma, St. Louis, MO, USA). The mice were euthanized 8 hr after hCG treatment and ovaries harvested.

Staging of estrous cycle by vaginal cytology

Request a detailed protocol

As previously described (Kano et al., 2017; Byers et al., 2012), staging of mice was performed using a wet cotton swab, introduced into the vaginal orifice then smeared onto a glass slide which was air-dried, stained with Giemsa, and scored for cytology by two independent observers. Briefly, proestrus was determined if the smear showed a preponderance of nucleated epithelial cells as well as leukocytes. Estrous was marked by an abundance of cornified epithelial cells, while metestrous smears contained a mixture of cornified epithelial cells and leukocytes. Finally, diestrus was characterized by abundant leukocytes with low numbers of cornified epithelium or nucleated epithelial cells.

Generation of single-cell suspension

Request a detailed protocol

Single-cell suspension from mouse ovaries was obtained as previously described with uterine enzymatic dissociation (Saatcioglu et al., 2019). Briefly, ovaries were incubated for 30 min at 34°C in dissociation medium (82 mM Na2SO4, 30 mM K2SO4, 10 mM Glucose, 10 mM HEPES, and 5 mM MgCL2, pH 7.4) containing 15 mg of Protease XXIII (Worthington), 100 U Papain, with 5 mm L-Cysteine, 2.5 mM EDTA (Worthington), and 1333 U of DNase 1 (Worthington). The reaction was then stopped in cold medium, and samples were mechanically dissociated, filtrated, and spun down three times before being resuspended to a concentration of 150,000 cells/mL in 20% Optiprep (Sigma) for inDrop sorting.

Single-cell RNA sequencing (inDrop)

Request a detailed protocol

Fluidic sorting was performed using the inDrop platform at the Single-Cell Core facility at Harvard Medical School as previously described (Klein et al., 2015; Macosko et al., 2015). We generated libraries of approximately 1500 cells per animal which were sequenced using the NextSeq500 (Illumina) platform. Transcripts were processed according to a previously published pipeline Klein et al., 2015 used to build a custom transcriptome from the Ensemble GRCm38 genome and GRCm38.84 annotation using Bowtie 1.1.1. Unique molecular identifiers (UMIs) were used to reference sequence reads back to individual captured molecules, referred to as UMIFM counts. All steps of the pipeline were run using default parameters unless explicitly specified.

scRNAseq data analysis

Data processing

Request a detailed protocol

The initial Seurat object was created using thresholds to identify putative cells (unique cell barcodes) with the following parameters: 1000–20,000 UMIs, 500–5000 genes, and less than 15% mitochondrial genes. The final merged dataset contained ~70,000 cells which were clustered based on expression of marker genes. These were further processed in several ways to exclude low-quality data and potential doublets. Visualization of single-cell data was performed using a non-linear dimensionality-reduction technique, uniform manifold approximation and projection. Markers for each level of cluster were identified using MAST in Seurat (R version 4.1.3 - Seurat version 4.1.0). Following identification of the main clusters (granulosa, mesenchyme, endothelium, immune, epithelium, and oocyte), we reanalyzed each cluster population to perform subclustering. Briefly, the granulosa, mesenchyme, and epithelium clusters were extracted from the integrated dataset by the subset function. The isolated cluster was then divided into several subclusters following normalization, scale, principal component analysis (PCA), and dimensionality reductions as previously described (Niu and Spradling, 2020).

Volcano plots

Request a detailed protocol

Highly differentially expressed genes between different estrous cycles were identified using the function FindMarkers in Seurat. Volcano plots were generated using ggplot2 package in R.

Pathway enrichment analysis

Request a detailed protocol

Differentially expressed genes with at least twofold changes between contiguous estrous stages were used as input for gene ontology enrichment analysis by clusterProfiler. Enrichplot package was used for visualization. Biological process subontology was chosen for this analysis.

Principal component analysis

Request a detailed protocol

PCA was used to identify common patterns of gene expression across stages of the cycle. For each Level 0 cluster object, cycling cells were extracted, and genes that were expressed in more than 5% of cells were identified. The expression of these genes in the cycling cells were scaled (set to mean zero, SD 1) and averaged across each of the four cycle stages. PCA was run (prcomp) on the average scaled expression data.

In situ hybridization and immunohistochemistry

Request a detailed protocol

In situ hybridizations were performed using ACDBio kits as per manufacturer’s protocol, as previously described (Saatcioglu et al., 2019). Briefly, RNAish was developed using the RNAscope 2.5 HD Reagent Kit (RED and Duplex, ACD Bio). Following deparaffinization in xylene, dehydration, peroxidase blocking, and heat-induced epitope retrieval by the target retrieval and protease plus reagents (ACD bio), tissue sections were hybridized with probes for the target genes (see Key resources table for accession number and catalog number of each gene) in the HybEZ hybridization oven (ACD Bio) for 2 hr at 40°C. The slides were then processed for standard signal amplification steps and chromogen development. Slides were counterstained in 50% hematoxylin (Dako), air dried, and cover-slipped with EcoMount. In addition to cycling and non-cycling mice, superovulated mice were used to validate markers from follicles associated with LH surge response in ovulatory follicles at the estrous stage for more precise timing of collection.

For colocalization of RNAish staining with immunohistochemistry, we first processed the tissue section for RNAscope as described above, including deparaffinization, antigen retrieval, hybridization, and chromogen development. Sections were then blocked in 3% bovine serum albumin in Tris-buffered solution (TBS) for 1 hr. Following three washes with TBS, the sections were incubated with the primary antibody (smooth muscle actin primary antibody; 1:300, Abcam) overnight at 4°C and developed with Dako EnVision + System horseradish peroxidase (HRP). Labeled polymer anti-rabbit was used as the secondary antibody, and the HRP signal was detected using the Dako detection system. Slides were then counterstained in hematoxylin and mounted as described above.

Reverse transcription-quantitative polymerase chain reaction

Request a detailed protocol

Mice were monitored through the estrous cycle and sacrificed at specific stage/timepoints as described above. Ovaries were dissected, and total RNA was extracted using the Qiagen RNA extraction kit (Qiagen). A cDNA library was synthesized from 500 ng total RNA using SuperScript III First-Strand Synthesis System for RT-PCR using manufacturer’s instructions with random hexamers (Invitrogen). The primers used for this study are described in Supplementary file 1. Expression levels were normalized to the Gapdh transcript using cycle threshold (Ct) values logarithmically transformed by the 2−ΔCt function.

ELISA

Request a detailed protocol

Blood was collected from mice by facial vein puncture, incubated at room temperature (RT) until spontaneously clotted, centrifuged at 8000 rpm for 5 min to collect the serum layer, and diluted 1/10 in each ELISA kit according to the manufacturing protocol; Mouse CNP/NPPC ELISA kit; Mouse serine protease inactive 35 (PRSS35) ELISA kit; Mouse TINAGL1 /Lipocalin 7 ELISA kit; and Human/Mouse/Rat Activin A Quantikine ELISA Kit (see Key resources table).

Data availability

Sequencing data have been deposited in Open Science Framework and the Broad Institute Single Cell Portal under study number SCP1914.

The following data sets were generated
    1. David P
    (2022) Open Science Framework
    ID 924fz. Single Cell Sequencing of the Mouse Ovary in Diverse Reproductive States.

References

    1. Dunlop CE
    2. Anderson RA
    (2014) The regulation and assessment of follicular growth
    Scandinavian Journal of Clinical and Laboratory Investigation. Supplementum 244:13–17.
    https://doi.org/10.3109/00365513.2014.936674

Decision letter

  1. Valerie Horsley
    Reviewing Editor; Yale University, United States
  2. Marianne E Bronner
    Senior Editor; California Institute of Technology, United States

Our editorial process produces two outputs: (i) public reviews designed to be posted alongside the preprint for the benefit of readers; (ii) feedback on the manuscript for the authors, including requests for revisions, shown below. We also include an acceptance summary that explains what the editors found interesting or important about the work.

Decision letter after peer review:

Thank you for submitting your article "A single cell atlas of the cycling murine ovary" for consideration by eLife. Your article has been reviewed by 3 peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Marianne Bronner as the Senior Editor. The reviewers have opted to remain anonymous.

The reviewers have discussed their reviews with one another, and the Reviewing Editor has drafted this to help you prepare a revised submission.

Essential revisions:

The major concerns were regarding some data analysis and a request for further validation of the RNA seq. data.

1. In the Monocle in Figure 5C-D, why are the periovulatory GC that is the visualization look very much like mural GC in a completely different branch than mural GC? I have the impression that either the clusters need to be refined/validated further or the Monocle is not providing meaningful results regarding lineage trajectories. Why would the "periovulatory" GCs be at the end of the CL branch? In addition, the population of periovulatory Runx1+ cells seems rather aspecific. I think the authors could increase the robustness of the data/evidence regarding this cluster. Particularly because they go on comparing periovulatory and antral GCs. Could at least one differentially expressed marker from the volcano plot be validated?

2. The author used tSNE for visualization of the generated scRNA seq dataset, which, according to my knowledge, is outdated for scRNA seq data visualization as its reproducibility has become an issue. Which version of the Seurat package does the author use? And also the other software information should be implemented. Therefore, I suggest the author reanalyze their dataset using an updated Seurat pipeline, and also reanalyze all of their data using UMAP.

3. It does not appear that Figure 2 represents a real "subcluster" analysis. Instead, these cells are from the original clustering performed on all cells. To get a better resolution as to the different cell types, it would be best to perform a subcluster analysis where these cells are first extracted from the total data set and then re-clustered using more relevant principal components (for example, see Niu and Spradling, 2020).

4. For the pseudotime analysis presented in Figure 4 the authors excluded the mitotic and putative atretic granulosa cells. What is the justification for this?

5. The reviewers also have many other suggestions that the authors should likely consider to clarify the manuscript for the broad audience of eLife readers.

Reviewer #1 (Recommendations for the authors):

I think this is a fantastic dataset, but I think it could be better explored. The relation between different stages of estrous/lactating and the histology/validation could have been better explored and the cell types more explicitly named. I am not sure the Monocle analysis is clarifying things, this is tricky to use and I would consider using other comparable methods to determine the trajectories.

I have some points that I would like to see clarified:

– How was it determined that the postpartum 10-day females were lactating or not-lactating?

– Could the authors color the cells per oestrus phase, postlactating, randomly next to Figure 1B? It would be very useful to understand how the different cells in the initial clustering are represented in the different groups, particularly to determine the periovulatory follicles. In that sense, it would be very helpful to have corresponding histology of the ovaries in the respective cycle phase from the beginning. This aspect of the manuscript should be better explored. Also, the ovary from the lactating/nonlactating could be provided.

– Could the authors clarify what cluster corresponds to theca externa? In the text, they mentioned the expression of Mfap5, but in Figure 2 they only mention smooth muscle cells. How are the population for the smooth muscle cells and theca externa connected and or differentiated? This is confusing the text and Figures. In other words: where are the theca externa cells in Figure 2? Perhaps the authors could provide fields with higher magnification, eventually in the Suppl Figure? Could the "smooth muscle cells" simply be the theca externa?

– In Figure 2 is also not clear where are the theca interna cells, although the authors do mention those in the text. This is confusing. Please refer to the theca interna in the text and Figure. Could there be two populations of Theca interna (steroidogenic) and the ones closer to the basement membrane? I don't see any evidence to call this 'immature' theca (see Figure 4C).

– Regarding the atretic granulosa cells: the follicle provided still looks rather healthy. Is there a second marker you could use together with Ghr to confirm that particular cluster corresponds to atretic GC? Or can you colocalise Ghr with a marker of atresia?

– In addition, the periovulatory marker chosen does not seem very specific. I would suggest the authors choose a better marker. Otherwise, the identity of this cluster remains inconclusive. Moreover, what is the evidence that this would be periovulatory follicles? Is this the phase of the cycle? This could be better explained/evidenced, perhaps using histology of ovaries corresponding to the respective stages.

– The 3 clusters of granulosa cells (Figure 3) in the corpus luteum are intriguing. I wonder whether the authors have an explanation for that. What are the oestrous stages the cells were observed in?

– I am not sure the Monocle results (Figure 4A-B) are informative because the theca are not separated in theca interna and externa. Hence to me, it is unclear how the mesenchymal clusters, if they are well determined and reflect biology, relate to each other.

– In the Monocle in Figure 5C-D, I have similar issues – why are the periovulatory GC that is the visualization look very much like mural GC in a completely different branch than mural GC? I have the impression that either the clusters need to be refined/validated further or the Monocle is not providing meaningful results regarding lineage trajectories. Why would the "periovulatory" GCs be at the end of the CL-branch? In addition, the population of periovulatory Runx1+ cells seems rather aspecific. I think the authors could increase the robustness of the data/evidence regarding this cluster. Particularly because they go on comparing periovulatory and antral GCs. Could at least one differentially expressed marker from the volcano plot be validated?

– I am not sure the Monocle analysis is clarifying things, this is tricky to use and I would consider using other comparable methods to determine the trajectories.

– In Figure 4C: if the authors claim that the stroma is more medullary or cortex. In the Results, they do not visualise it. Hence, I would remove the claim from the Results but keep it in the Discussion.

– Results on the OSE should be shown in the Results section, now they are part of the Discussion? Where are those cells? What clusters? Etc. This population is very interesting/elusive and should be given proper attention/validation in the Results. In addition, the authors mention no differences between the estrous states: this should be clearly shown in the Results.

– The Discussion is too long, in particular, the discussion on the biomarkers reads more like a review than a discussion of your results. This should be shortened and made concise.

Reviewer #3 (Recommendations for the authors):

1. The introduction is very brief. It would be improved with a more thorough description of the estrus cycle beyond simply naming the stages. For example, how many days is the estrus cycle in the mouse? This would make the paper more accessible to non-experts.

2. Data availability- The authors have deposited the count matrix of the combined inDROP scRNA-seq dataset at Open Science Framework. However, for this dataset to serve as a useful resource for the wider community, the processed data should be uploaded to the Broad Institute Single Cell Portal (https://singlecell.broadinstitute.org).

3. The abstract indicates that in addition to ovaries isolated from the four stages of estrus, scRNA-seq was also performed on ovaries isolated from lactating or non-lactating 10 days postpartum mice, and from randomly cycling mice. However, cells from these later samples were never presented or discussed. It is therefore not clear what cells were used in the analysis. To make this clear, Figure 1 (or a supplemental figure) should include a tSNE plot where the cells are color-coded based on their sample origin (i.e. like Figure 5A).

4. Figure 1 – Given that this is the first figure, it would be good to show a tSNE plot that identifies which cells are from which library.

5. Figure 2 – It does not appear that Figure 2 represents a real "subcluster" analysis. Instead, these cells are from the original clustering performed on all cells. To get a better resolution as to the different cell types, it would be best to perform a subcluster analysis where these cells are first extracted from the total data set and then re-clustered using more relevant principal components (for example, see Niu and Spradling, 2020).

6. For the data presented in Figure 2, it is not entirely clear what is novel and what is already known. There are no references to the genes used to identify the cell sub-types types. References that are cited in the tables are not listed in the reference list. These references should, at a minimum, be listed in a supplemental text and it should be made clear in the text what genes are novel.

7. In addition to the markers chosen for specificity (e.g. Mafap5), it would be nice to see representative UMAPS (or dot plots) showing standard cell-type markers (e.g. Tcf21, Dcn, Notch3, Cxcl12). This could be a supplemental figure.

8. The description of the mesenchyme cell subpopulations is confusing as the authors do not use consistent terminology. A more concise definition of theca interna, theca externa, mature theca and immature theca as it relates to theca 1 and 2, stroma 1 and 2, smooth muscle would make the paper more accessible to the non-expert.

9. Line 196: Text refers to "dividing mesenchyme (8%) as seen in Figure 2A" but this population is not labeled. Which subcluster of cells is this referring to?

10. Line 204: Confusing to use theca externa when Figure 2C is labeled smooth muscle. Defining theca interna vs. theca externa, as mentioned above, would help.

11. Line 206: "…whereas Hhip was expressed in theca interna and immature theca (Figure S2A, B)." What is the difference between theca 1 and immature theca? Are these the same?

12. Figure 2 FigSup2: hard to see Acta2-Mfap5 overlap in Figure S2B. A higher magnification image would be helpful. It looks like there are two distinct cell layers, not the same cells. Also, Hhip expressing cells appear to be located outside of the Acta2-expressing cells, but Hhip is referred to as marking theca interna while Acta2 is said to mark theca externa. This needs to be explained.

13. Figure 2C – To make interpretation of the expression patterns more accessible to the non-expert, I suggest labeling some structures in these figure panels (e.g. oocyte, granulosa cells, cortical, medullary).

14. What cell type are the mitotic granulosa cells most similar to? Antral/mural or preantral/cumulus? Or are they an intermediate between these two? Subcluster analysis, as mentioned above, may help to better define their relationship.

15. Line 228: What is the justification for identifying these as atretic? What is the significance of Ghr to atresia? References would help.

16. Fan et al., used the following gene expression signatures to distinguish cumulus from mural GC's: cumulus GC (VCANhigh/FSThigh/IGFBP2high/HTRA1high/INHBBhigh/IHHhigh); mural GC (WT1low/EGR4low/KRT18high/CITED2high/LIHPhigh/AKIRIN1high). However, in Table S5 Fst and Inhbb are listed as markers for mural, not cumulus cells. This should be explained.

17. For the pseudotime analysis presented in Figure 4 the authors excluded the mitotic and putative atretic granulosa cells. What is the justification for this?

18. Figure 4D and E – Why are the preovulatory cells at the terminus of the pseudotime with CL2, while CL1 and CL3 cells positioned along the root. This does not match expectations. The discussion suggests that the developmental progression is PO>CL2>CL1>CL3, but this is not supported by the trajectory analysis. This calls into question the usefulness of the pseudotime analysis. It is mentioned in the discussion that instead of the different CL clusters representing a developmental progression, they are instead distinct cell types within a CL. This could be resolved with double RNA in situ hybridization using CL cluster-specific genes.

19. Figure 5A – It would be nice to show 4 separate plots for each of the stages. Hard to see this on a single plot. Perhaps 4 smaller panels next to A. It would also be helpful to label the different cell sub-types so the reader does not need to refer back to Figure 3.

20. Figure 5B – It would be helpful to label this panel "Proestrus-Estrus," following that in Figure S5.

21. Line 270 – reference to Figure S5D should be Figures 5E.

22. Figure S5D – This figure panel is not mentioned in the text.

23. Line 270: "…involved in wound repair during ovulation." This needs a reference.

24. It is stated that the motivation for identifying stage-specific secreted biomarkers was to identify markers that would be useful for staging in assisted reproduction and other applications in reproductive medicine. This begs the question if Prss35, Nppc and Tinagl1 are also expressed in human ovaries? Was this part of the reason they were selected?

25. Figure 6B: I think panel B should be first since it is from the scRNA-seq data set and then A is validation. I would also add Lhcgr and Pgr to the dot plot in B. Also, it makes more sense to put Inhba next to Nppc.

26. Figure 6A: The graphs in A are somewhat randomly organized, which makes it unnecessarily complicated. These should be reorganized to group factors high in P vs. E together and those high in E vs. P together.

27. Figure 6C: Follicle stages should be labeled. Mural vs. Cumulus cell should be labeled. Would be best to show both follicles that are labeled and ones that are not labeled to emphasize that this is a stage-specific expression. I cannot see Nppc expression. Might be helpful to have a higher magnification or arrows pointing to expressing cells.

28. Line 234: Run-on sentence: "Early pre-antral follicles…"

29. Line 384- Fan et al., reference should be deleted here.

https://doi.org/10.7554/eLife.77239.sa1

Author response

Essential revisions:

The major concerns were regarding some data analysis and a request for further validation of the RNA seq. data.

1. In the Monocle in Figure 5C-D, why are the periovulatory GC that is the visualization look very much like mural GC in a completely different branch than mural GC? I have the impression that either the clusters need to be refined/validated further or the Monocle is not providing meaningful results regarding lineage trajectories. Why would the "periovulatory" GCs be at the end of the CL branch?

We agree that the pseudo-time analysis we presented using monocle appeared counter intuitive. Furthermore, other pseudotime analyses gave similar results. This is because the lineage from preantral granulosa cell to cumulus granulosa cells of antral follicles is linear. In contrast, mural granulosa cells of antral follicles have a distinct trajectory representing a branching point leading to luteal differentiation. However to luteinizing cell state is unique, and more similar to luteal cells of the corpus luteum than to antral mural granulosa. This creates a discontinuous trajectory of granulosa lineages in pseudotime, where the transcriptional response to the LH surge appears as a terminal state, when it is in fact transient. Therefore, in an effort to improve clarity, we have removed the pseudotime analysis in favor of more concise explanation of these lineages and states in the discussion.

In addition, the population of periovulatory Runx1+ cells seems rather aspecific. I think the authors could increase the robustness of the data/evidence regarding this cluster. Particularly because they go on comparing periovulatory and antral GCs.

To improve the robustness of these conclusions we have replaced Runx1 by Oxtr which is more specific to this cell type (See Figure 3D and 3-supplement 1D). In addition, we refer to a swath of previously validated markers which are very specific to Graafian follicles.

Could at least one differentially expressed marker from the volcano plot be validated?

We validated several of the markers shown in the volcano plot as presented in in Fig5E.

2. The author used tSNE for visualization of the generated scRNA seq dataset, which, according to my knowledge, is outdated for scRNA seq data visualization as its reproducibility has become an issue. Which version of the Seurat package does the author use? And also the other software information should be implemented. Therefore, I suggest the author reanalyze their dataset using an updated Seurat pipeline, and also reanalyze all of their data using UMAP.

Following the reviewers’ advice, we reanalyzed our dataset and presented the dimensional reductions using UMAP. Information regarding the Seurat package and R version have been added to the material and methods section.

3. It does not appear that Figure 2 represents a real "subcluster" analysis. Instead, these cells are from the original clustering performed on all cells. To get a better resolution as to the different cell types, it would be best to perform a subcluster analysis where these cells are first extracted from the total data set and then re-clustered using more relevant principal components (for example, see Niu and Spradling, 2020).

When reanalyzing all our data with UMAP we re-clustered each cell type independently, as Niu and Spradling.

4. For the pseudotime analysis presented in Figure 4 the authors excluded the mitotic and putative atretic granulosa cells. What is the justification for this?

Given the reviewers’ comments we have elected to remove the pseudotime trajectory analysis in favor of more concise interpretation of lineages.

5. The reviewers also have many other suggestions that the authors should likely consider to clarify the manuscript for the broad audience of eLife readers.

We thank the reviewers for their constructive comments regarding the manuscript. We have edited the presentation of our findings to make it more accessible to the broad readership of eLife and anticipate that this manuscript will provide a useful resource for the research community regarding the dynamic changes in transcriptome at the single cell level in the ovary as a function of the estrous cycle.

Reviewer #1 (Recommendations for the authors):

I think this is a fantastic dataset, but I think it could be better explored. The relation between different stages of estrous/lactating and the histology/validation could have been better explored and the cell types more explicitly named. I am not sure the Monocle analysis is clarifying things, this is tricky to use and I would consider using other comparable methods to determine the trajectories.

I have some points that I would like to see clarified:

– How was it determined that the postpartum 10-day females were lactating or not-lactating?

The pups were removed from the female on the day they were born, hence when the ovaries were harvested at post-partum day 10 while the mice were not-lactating. The lactating mice remained with their pups at post-partum day 10. The text has been modified to clarify the method (lines 77-80).

– Could the authors color the cells per oestrus phase, postlactating, randomly next to Figure 1B? It would be very useful to understand how the different cells in the initial clustering are represented in the different groups, particularly to determine the periovulatory follicles. In that sense, it would be very helpful to have corresponding histology of the ovaries in the respective cycle phase from the beginning. This aspect of the manuscript should be better explored. Also, the ovary from the lactating/nonlactating could be provided.

We added a dimensional reduction showing the distribution of the cells depending on the stage of the estrous cycle in which they have been collected (Figure 1B). The corresponding representative histology for each phase of the estrous cycle as well as day 10 post-partum lactating and non-lactating is illustrated on Figure S1B.

– Could the authors clarify what cluster corresponds to theca externa? In the text, they mentioned the expression of Mfap5, but in Figure 2 they only mention smooth muscle cells. How are the population for the smooth muscle cells and theca externa connected and or differentiated? This is confusing the text and Figures. In other words: where are the theca externa cells in Figure 2? Perhaps the authors could provide fields with higher magnification, eventually in the Suppl Figure? Could the "smooth muscle cells" simply be the theca externa?

We added some clarification of the cell types present in the theca externa which includes the Smooth muscle cells which we validated in greater detail (lines 220-229). We also provided in supplemental figure 2 – supplement 1 a higher magnification of Hhip-Acta2 and well as Mfap5-Acta2 to be able to differentiate smooth muscle cells of the theca externa from the early theca cells of the theca interna.

– In Figure 2 is also not clear where are the theca interna cells, although the authors do mention those in the text. This is confusing. Please refer to the theca interna in the text and Figure. Could there be two populations of Theca interna (steroidogenic) and the ones closer to the basement membrane? I don't see any evidence to call this 'immature' theca (see Figure 4C).

The theca cell clusters were renamed Steroidogenic and early theca. We provided a new micrograph of Hhip in situ hybridization to show expression in “early” theca cells surrounding preantral follicles, while Cyp17a1 is expressed in “steroidogenic” theca cells surrounding antral follicles in Figures 2C and Figure 2 – supplement 1A, B.

– Regarding the atretic granulosa cells: the follicle provided still looks rather healthy. Is there a second marker you could use together with Ghr to confirm that particular cluster corresponds to atretic GC? Or can you colocalise Ghr with a marker of atresia?

We provided a better representative example of atretic follicle. To confirm that this cluster corresponds to atretic granulosa cells, we provided in Figure S3A a feature plots of other markers known to be expressed in apoptotic cells that colocalized with Ghr.

– In addition, the periovulatory marker chosen does not seem very specific. I would suggest the authors choose a better marker. Otherwise, the identity of this cluster remains inconclusive. Moreover, what is the evidence that this would be periovulatory follicles? Is this the phase of the cycle? This could be better explained/evidenced, perhaps using histology of ovaries corresponding to the respective stages.

We replaced Runx1 by Oxtr as the representative periovulatory marker, which appears more specific by in-situ hybridization (Figures3C and S3D). The classification of this cell type is further supported by a fullyreferenced list of markers in Table S5. Furthermore, this cluster is only found in mice as estrous as indicated in Figure 5A/B, which is consistent spatially and temporally with mural granulosa cells of periovulatory graafian follicles responding to LH and initiating luteinization.

– The 3 clusters of granulosa cells (Figure 3) in the corpus luteum are intriguing. I wonder whether the authors have an explanation for that. What are the oestrous stages the cells were observed in?

Following reanalysis of our dataset with the latest version of Seurat and the dimensional representation method to UMAP we now have two corpus luteum (CL) clusters instead of 3. We hypothesize that these two clusters correspond to active and regressing CL respectively. We confirmed this hypothesis by evaluating the level of expression of a proliferative marker such as Top2a and a cell cycle arrest such as Cdkn1a (Figure S3B). We further validated the spatial colocalization of these markers, with large healthy corpora lutea expressing Top2a while smaller luteolysing CLs express the cell cycle inhibitor and senescence marker Cdkn1a along with a cohort of validated markers shown in FigS3E and Table S4.

– I am not sure the Monocle results (Figure 4A-B) are informative because the theca are not separated in theca interna and externa. Hence to me, it is unclear how the mesenchymal clusters, if they are well determined and reflect biology, relate to each other.

We have removed the pseudotime analysis from the manuscript.

– In the Monocle in Figure 5C-D, I have similar issues – why are the periovulatory GC that is the visualization look very much like mural GC in a completely different branch than mural GC? I have the impression that either the clusters need to be refined/validated further or the Monocle is not providing meaningful results regarding lineage trajectories.

We agree with the reviewer. This figure was removed.

Why would the "periovulatory" GCs be at the end of the CL-branch? In addition, the population of periovulatory Runx1+ cells seems rather aspecific. I think the authors could increase the robustness of the data/evidence regarding this cluster. Particularly because they go on comparing periovulatory and antral GCs. Could at least one differentially expressed marker from the volcano plot be validated?

As discussed above, the position of the periovulatory granulosa cells in the monocle visual representation of pseudotime is a function of the LH response signature state, rather than a linear differentiation trajectory, which complicates the interpretation of those plots, as the cells then revert along the same axis to a luteinized state with reduced acute response to LH. To clarify the manuscript, we have removed the pseudotime plot. The classification of periovulatory GC is further supported by a fully-referenced list of markers in Table S5.

– I am not sure the Monocle analysis is clarifying things, this is tricky to use and I would consider using other comparable methods to determine the trajectories.

We agree with the reviewer. This figure was removed.

– In Figure 4C: if the authors claim that the stroma is more medullary or cortex. In the Results, they do not visualise it. Hence, I would remove the claim from the Results but keep it in the Discussion.

We agree with the reviewer. Following reanalyzing our data with UMAP we couldn’t identify unique dichotomous markers to differentiate each cluster. However, we found Ennp2 to be differentially expressed: present in one cluster (fibroblast-like stroma) and completely absent from the other (steroidogenic stroma). These results are now part of figure 2 and figure 2 – supplement 1.

– Results on the OSE should be shown in the Results section, now they are part of the Discussion? Where are those cells? What clusters? Etc. This population is very interesting/elusive and should be given proper attention/validation in the Results. In addition, the authors mention no differences between the estrous states: this should be clearly shown in the Results.

The epithelium cluster is composed of the ovarian surface epithelium cells. We now explore this cluster deeper in Figure 4. We showed that this cluster could be subdivided in two clusters one of them being characterized by a very strong expression of proliferation markers and being composed almost exclusively of granulosa cells collected during the estrous phase of the cycle, consistent with a transient amplification of these cells coinciding with ovulatory wound repair.

– The Discussion is too long, in particular, the discussion on the biomarkers reads more like a review than a discussion of your results. This should be shortened and made concise.

We thank the reviewer for this suggestion, we significantly revised the discussion to clarify our results in the context of ovarian biology and to make it more concise.

Reviewer #3 (Recommendations for the authors):

1. The introduction is very brief. It would be improved with a more thorough description of the estrus cycle beyond simply naming the stages. For example, how many days is the estrus cycle in the mouse? This would make the paper more accessible to non-experts.

We thank the reviewer for this suggestion, we significantly rewrote the introduction and discussion to improve clarity and make the manuscript more accessible to the broad readership of eLife. We hope this dataset will become a valuable resource to the community to investigate dynamic cell states associated with estrous cycling in the ovary.

2. Data availability- The authors have deposited the count matrix of the combined inDROP scRNA-seq dataset at Open Science Framework. However, for this dataset to serve as a useful resource for the wider community, the processed data should be uploaded to the Broad Institute Single Cell Portal (https://singlecell.broadinstitute.org).

We have deposited the data set and source code on the Broad Institute Single Cell Portal under study number SCP1914.

3. The abstract indicates that in addition to ovaries isolated from the four stages of estrus, scRNA-seq was also performed on ovaries isolated from lactating or non-lactating 10 days postpartum mice, and from randomly cycling mice. However, cells from these later samples were never presented or discussed. It is therefore not clear what cells were used in the analysis. To make this clear, Figure 1 (or a supplemental figure) should include a tSNE plot where the cells are color-coded based on their sample origin (i.e. like Figure 5A).

We added panels to Figure 1B and 5A/B to represent the distribution of cells depending on the stage of the estrous cycle. We also illustrated the morphology of the ovary depending on the stage of the estrous cycle and lactation state on Figure S1B. Finally, we analyzed the post-partum lactation data in Figure S3E and discuss it in the results (lines 272-275).

4. Figure 1 – Given that this is the first figure, it would be good to show a tSNE plot that identifies which cells are from which library.

We now reran the analysis to show UMAP dimensional reduction instead of tSNE and updated all figures of the manuscript. We also added a panel in Figure 1B to identify which cell are coming from which library.

5. Figure 2 – It does not appear that Figure 2 represents a real "subcluster" analysis. Instead, these cells are from the original clustering performed on all cells. To get a better resolution as to the different cell types, it would be best to perform a subcluster analysis where these cells are first extracted from the total data set and then re-clustered using more relevant principal components (for example, see Niu and Spradling, 2020).

We thank the reviewer for this suggestion. We re-clustered each object and displayed the results using UMAP projections and revised all the figures of the manuscript accordingly.

6. For the data presented in Figure 2, it is not entirely clear what is novel and what is already known. There are no references to the genes used to identify the cell sub-types types. References that are cited in the tables are not listed in the reference list. These references should, at a minimum, be listed in a supplemental text and it should be made clear in the text what genes are novel.

We thank the reviewer for pointing this out. The top markers for each cluster/cell type are presented in Tables S2, S3, and S4. The markers that have previously been described are fully referenced in these tables, and those citations are now part of the manuscript’s bibliography.

7. In addition to the markers chosen for specificity (e.g. Mafap5), it would be nice to see representative UMAPS (or dot plots) showing standard cell-type markers (e.g. Tcf21, Dcn, Notch3, Cxcl12). This could be a supplemental figure.

We thank the reviewer for this suggestion. A DotPlot was used to show the expression of the suggested and validated cell-type markers, which is now located in Figure S2D.

8. The description of the mesenchyme cell subpopulations is confusing as the authors do not use consistent terminology. A more concise definition of theca interna, theca externa, mature theca and immature theca as it relates to theca 1 and 2, stroma 1 and 2, smooth muscle would make the paper more accessible to the non-expert.

The new clustering analysis and nomenclature clarify the mesenchymal cell types. The description of the theca interna and externa, and their cellular composition is now discussed more clearly in the results (lines 230-244) and discussion (lines 297- 417), and visualized in Figure 1 – supplement 1. We also added a zoomed picture of the theca stained with SMA and Mfap5 or Hhip to help the reader differentiate smooth muscles of the theca externa to the early theca cells of the theca interna depending on expression and position in figure 2 – supplement 1B.

9. Line 196: Text refers to "dividing mesenchyme (8%) as seen in Figure 2A" but this population is not labeled. Which subcluster of cells is this referring to?

This population no longer represent a distinct cluster in our new analysis.

10. Line 204: Confusing to use theca externa when Figure 2C is labeled smooth muscle. Defining theca interna vs. theca externa, as mentioned above, would help.

The description of the theca interna and externa, and their cellular composition is now discussed more clearly in the results (lines 230-244) and discussion (lines 397-417), and visualized in Figure 1supplement 1A

11. Line 206: "…whereas Hhip was expressed in theca interna and immature theca (Figure S2A, B)." What is the difference between theca 1 and immature theca? Are these the same?

We adopted a consistent nomenclature reflecting the updated subclustering analysis. The theca clusters have been renamed early theca, steroidogenic theca, and smooth muscle (of the theca externa) and the text was modified accordingly.

12. Figure 2 FigSup2: hard to see Acta2-Mfap5 overlap in Figure S2B. A higher magnification image would be helpful. It looks like there are two distinct cell layers, not the same cells. Also, Hhip expressing cells appear to be located outside of the Acta2-expressing cells, but Hhip is referred to as marking theca interna while Acta2 is said to mark theca externa. This needs to be explained.

We added higher magnification of Hhip and Mfap5 staining colocalization with Acta2 in Figure 2 – supplement 1B. While Hhip is expressed in the early theca, of the theca interna, Mfap5 is expressed in smooth muscle cells of the theca externa and colocalizes with Acta2, a general marker of smooth muscle.

13. Figure 2C – To make interpretation of the expression patterns more accessible to the non-expert, I suggest labeling some structures in these figure panels (e.g. oocyte, granulosa cells, cortical, medullary).

We thank the reviewer for this suggestion. The general histology of ovarian structures across diverse reproductive states are now shown in supplemental figure 1 – supplement 1A.

14. What cell type are the mitotic granulosa cells most similar to? Antral/mural or preantral/cumulus? Or are they an intermediate between these two? Subcluster analysis, as mentioned above, may help to better define their relationship.

The mitotic granulosa cells cluster was composed of cells that were closely related to both the antral mural and the preantral-cumulus granulosa (see Author response image 2 the expression of both markers Mro (red) and Kctd14(green) in the mitotic cluster). While these clusters may have been differentiated with a higher cluster resolution, we chose to keep the resolution consistent across all subclustering.

Author response image 1

15. Line 228: What is the justification for identifying these as atretic? What is the significance of Ghr to atresia? References would help.

We propose Ghr as a new marker for atretic granulosa cells. We confirmed the identification of atretic cluster based on expression of apoptotic markers such as Pik3ip1, Nupr1 and Gadd45a presented on Figure S3A along with additional previously validated markers which are referenced in Tables S4. Moreover, we provided a better representative example of an atretic follicle with Ghr in-situ hybridization in Figure 3C.

16. Fan et al., used the following gene expression signatures to distinguish cumulus from mural GC's: cumulus GC (VCANhigh/FSThigh/IGFBP2high/HTRA1high/INHBBhigh/IHHhigh); mural GC (WT1low/EGR4low/KRT18high/CITED2high/LIHPhigh/AKIRIN1high). However, in Table S5 Fst and Inhbb are listed as markers for mural, not cumulus cells. This should be explained.

While we investigated all markers discussed in the Fan et al., publication, many were not conserved from human to mouse (see Author response image 2). We elaborate on some of these inter-species’ discrepancies in the discussion (lines 414-417).

Author response image 2

17. For the pseudotime analysis presented in Figure 4 the authors excluded the mitotic and putative atretic granulosa cells. What is the justification for this?

Given that the monocle pseudotime trajectories of granulosa cells were difficult to interpret, we have removed this analysis from the manuscript.

18. Figure 4D and E – Why are the preovulatory cells at the terminus of the pseudotime with CL2, while CL1 and CL3 cells positioned along the root. This does not match expectations. The discussion suggests that the developmental progression is PO>CL2>CL1>CL3, but this is not supported by the trajectory analysis. This calls into question the usefulness of the pseudotime analysis. It is mentioned in the discussion that instead of the different CL clusters representing a developmental progression, they are instead distinct cell types within a CL. This could be resolved with double RNA in situ hybridization using CL cluster-specific genes.

We agree that the pseudo-time analysis we presented using monocle appears counter intuitive. This is because the pseudotime progression from mural granulosa cells of antral follicles to periovulatory follicle granulosa cells, to active CL to regressing CL does not appear linear. Therefore, in an effort to improve clarity, we have removed the pseudotime analysis in favor of more concise explanation of these lineages and states in the discussion (line 378-382).

19. Figure 5A – It would be nice to show 4 separate plots for each of the stages. Hard to see this on a single plot. Perhaps 4 smaller panels next to A. It would also be helpful to label the different cell sub-types so the reader does not need to refer back to Figure 3.

We thank the reviewer for this suggestion. We added a new figure (Figure 5B) showing 4 separate plots for each of the stages.

20. Figure 5B – It would be helpful to label this panel "Proestrus-Estrus," following that in Figure S5.

We thank the reviewer for bringing this to our attention, we have added a clear label.

21. Line 270 – reference to Figure S5D should be Figures 5E.

We thank the reviewer for bringing this to our attention, we fixed this mistake.

22. Figure S5D – This figure panel is not mentioned in the text.

We have added this to the text (line 317).

23. Line 270: "…involved in wound repair during ovulation." This needs a reference.

We have added references to the text (lines 317).

24. It is stated that the motivation for identifying stage-specific secreted biomarkers was to identify markers that would be useful for staging in assisted reproduction and other applications in reproductive medicine. This begs the question if Prss35, Nppc and Tinagl1 are also expressed in human ovaries? Was this part of the reason they were selected?

To prioritize these markers, we identified factors specifically secreted by ovarian granulosa in our dataset, and predominantly expressed by the gonads in the human GTEX database. We have added selection details to the text (lines 320-324).

25. Figure 6B: I think panel B should be first since it is from the scRNA-seq data set and then A is validation. I would also add Lhcgr and Pgr to the dot plot in B. Also, it makes more sense to put Inhba next to Nppc.

We thank the reviewer for the suggestion. We modified Figure 6B accordingly.

26. Figure 6A: The graphs in A are somewhat randomly organized, which makes it unnecessarily complicated. These should be reorganized to group factors high in P vs. E together and those high in E vs. P together.

We thank the reviewer for the suggestion. We modified Figure 6A accordingly.

27. Figure 6C: Follicle stages should be labeled. Mural vs. Cumulus cell should be labeled. Would be best to show both follicles that are labeled and ones that are not labeled to emphasize that this is a stage-specific expression. I cannot see Nppc expression. Might be helpful to have a higher magnification or arrows pointing to expressing cells.

We thank the reviewer for the suggestion. We labelled the structures in Fig6C and replaced the Nppc panel with a picture with a clearer stain.

28. Line 234: Run-on sentence: "Early pre-antral follicles…"

We thank the reviewer for bringing this to our attention, we fixed this mistake.

29. Line 384- Fan et al., reference should be deleted here.

We thank the reviewer for bringing this to our attention, we fixed this mistake.

https://doi.org/10.7554/eLife.77239.sa2

Article and author information

Author details

  1. Mary E Morris

    Department of Gynecology and Reproductive Biology, Massachusetts General Hospital, Boston, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Funding acquisition, Validation, Investigation, Writing - original draft
    Contributed equally with
    Marie-Charlotte Meinsohn
    For correspondence
    MEMORRIS@mgh.harvard.edu
    Competing interests
    No competing interests declared
  2. Marie-Charlotte Meinsohn

    1. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, Boston, United States
    2. Department of Surgery, Harvard Medical School, Boston, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Validation, Investigation, Methodology, Writing - original draft, Writing - review and editing
    Contributed equally with
    Mary E Morris
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0002-1378-4655
  3. Maeva Chauvin

    1. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, Boston, United States
    2. Department of Surgery, Harvard Medical School, Boston, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
  4. Hatice D Saatcioglu

    1. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, Boston, United States
    2. Department of Surgery, Harvard Medical School, Boston, United States
    Contribution
    Conceptualization, Data curation, Formal analysis, Investigation, Methodology
    Competing interests
    No competing interests declared
  5. Aki Kashiwagi

    1. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, Boston, United States
    2. Department of Surgery, Harvard Medical School, Boston, United States
    Contribution
    Data curation, Validation, Investigation, Methodology
    Competing interests
    No competing interests declared
  6. Natalie A Sicher

    1. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, Boston, United States
    2. Department of Surgery, Harvard Medical School, Boston, United States
    Contribution
    Data curation, Formal analysis, Validation, Investigation
    Competing interests
    No competing interests declared
  7. Ngoc Nguyen

    1. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, Boston, United States
    2. Department of Surgery, Harvard Medical School, Boston, United States
    Contribution
    Formal analysis, Validation, Investigation
    Competing interests
    No competing interests declared
  8. Selena Yuan

    1. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, Boston, United States
    2. Department of Surgery, Harvard Medical School, Boston, United States
    Contribution
    Data curation, Validation
    Competing interests
    No competing interests declared
  9. Rhian Stavely

    1. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, Boston, United States
    2. Department of Surgery, Harvard Medical School, Boston, United States
    Contribution
    Data curation, Formal analysis, Writing - review and editing
    Competing interests
    No competing interests declared
  10. Minsuk Hyun

    Howard Hughes Medical Institute, Department of Neurobiology, Harvard Medical School, Boston, United States
    Contribution
    Data curation, Software, Formal analysis, Investigation, Visualization
    Competing interests
    No competing interests declared
  11. Patricia K Donahoe

    1. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, Boston, United States
    2. Department of Surgery, Harvard Medical School, Boston, United States
    Contribution
    Supervision, Writing - original draft, Writing - review and editing
    Competing interests
    No competing interests declared
  12. Bernardo L Sabatini

    Howard Hughes Medical Institute, Department of Neurobiology, Harvard Medical School, Boston, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Investigation, Methodology, Writing - original draft, Project administration
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-0095-9177
  13. David Pépin

    1. Pediatric Surgical Research Laboratories, Massachusetts General Hospital, Boston, United States
    2. Department of Surgery, Harvard Medical School, Boston, United States
    Contribution
    Conceptualization, Data curation, Software, Formal analysis, Supervision, Funding acquisition, Investigation, Methodology, Writing - original draft, Project administration, Writing - review and editing
    For correspondence
    dpepin@mgh.harvard.edu
    Competing interests
    No competing interests declared
    ORCID icon "This ORCID iD identifies the author of this article:" 0000-0003-2046-6708

Funding

Eunice Kennedy Shriver National Institute of Child Health and Human Development (1R01HD102014-01)

  • David Pépin

Massachusetts General Hospital

  • David Pépin
  • Patricia K Donahoe

Huiying Foundation (Huiying Fellowship)

  • Hatice D Saatcioglu
  • David Pépin

The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.

Acknowledgements

We thank LiHua Zhang, Bugra Uluyurt, Phoebe A May, Caroline Coletti, and Sarah Mustafa Eisa for technical help. This study was supported by the National Institute for Child Health and Human Development to D.P. (1R01HD102014-01), the Huiying Fellowship (H.D.S. and D.P.), Sudna Gar Fellowship (D.P.), Massachusetts General Hospital Executive Committee on Research (D.P. and P.K.D.), and royalties (P.K.D.) from the use of the MIS ELISA in infertility clinics.

Ethics

Animal experiments were approved by the National Institute of Health and Harvard Medical School Institutional Animal Care and Use Committee, and performed in accordance with experimental protocols 2009N000033 and 2014N000275 approved by the Massachusetts General Hospital Institutional Animal Care and Use Committee.

Senior Editor

  1. Marianne E Bronner, California Institute of Technology, United States

Reviewing Editor

  1. Valerie Horsley, Yale University, United States

Publication history

  1. Received: January 20, 2022
  2. Preprint posted: February 11, 2022 (view preprint)
  3. Accepted: September 12, 2022
  4. Version of Record published: October 7, 2022 (version 1)

Copyright

© 2022, Morris, Meinsohn et al.

This article is distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use and redistribution provided that the original author and source are credited.

Metrics

  • 741
    Page views
  • 210
    Downloads
  • 0
    Citations

Article citation count generated by polling the highest count across the following sources: Crossref, PubMed Central, Scopus.

Download links

A two-part list of links to download the article, or parts of the article, in various formats.

Downloads (link to download the article as PDF)

Open citations (links to open the citations from this article in various online reference manager services)

Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)

  1. Mary E Morris
  2. Marie-Charlotte Meinsohn
  3. Maeva Chauvin
  4. Hatice D Saatcioglu
  5. Aki Kashiwagi
  6. Natalie A Sicher
  7. Ngoc Nguyen
  8. Selena Yuan
  9. Rhian Stavely
  10. Minsuk Hyun
  11. Patricia K Donahoe
  12. Bernardo L Sabatini
  13. David Pépin
(2022)
A single-cell atlas of the cycling murine ovary
eLife 11:e77239.
https://doi.org/10.7554/eLife.77239
  1. Further reading

Further reading

    1. Developmental Biology
    2. Genetics and Genomics
    Janani Ramachandran, Weiqiang Zhou ... Steven A Vokes
    Research Article Updated

    The larynx enables speech while regulating swallowing and respiration. Larynx function hinges on the laryngeal epithelium which originates as part of the anterior foregut and undergoes extensive remodeling to separate from the esophagus and form vocal folds that interface with the adjacent trachea. Here we find that sonic hedgehog (SHH) is essential for epithelial integrity in the mouse larynx as well as the anterior foregut. During larynx-esophageal separation, low Shh expression marks specific domains of actively remodeling epithelium that undergo an epithelial-to-mesenchymal transition (EMT) characterized by the induction of N-Cadherin and movement of cells out of the epithelial layer. Consistent with a role for SHH signaling in regulating this process, Shh mutants undergo an abnormal EMT throughout the anterior foregut and larynx, marked by a cadherin switch, movement out of the epithelial layer and cell death. Unexpectedly, Shh mutant epithelial cells are replaced by a new population of FOXA2-negative cells that likely derive from adjacent pouch tissues and form a rudimentary epithelium. These findings have important implications for interpreting the etiology of HH-dependent birth defects within the foregut. We propose that SHH signaling has a default role in maintaining epithelial identity throughout the anterior foregut and that regionalized reductions in SHH trigger epithelial remodeling.

    1. Developmental Biology
    Yanling Xin, Qinghai He ... Shuyi Chen
    Research Article

    N 6-methyladenosine (m6A) is the most prevalent mRNA internal modification and has been shown to regulate the development, physiology, and pathology of various tissues. However, the functions of the m6A epitranscriptome in the visual system remain unclear. In this study, using a retina-specific conditional knockout mouse model, we show that retinas deficient in Mettl3, the core component of the m6A methyltransferase complex, exhibit structural and functional abnormalities beginning at the end of retinogenesis. Immunohistological and single-cell RNA sequencing (scRNA-seq) analyses of retinogenesis processes reveal that retinal progenitor cells (RPCs) and Müller glial cells are the two cell types primarily affected by Mettl3 deficiency. Integrative analyses of scRNA-seq and MeRIP-seq data suggest that m6A fine-tunes the transcriptomic transition from RPCs to Müller cells by promoting the degradation of RPC transcripts, the disruption of which leads to abnormalities in late retinogenesis and likely compromises the glial functions of Müller cells. Overexpression of m6A-regulated RPC transcripts in late RPCs partially recapitulates the Mettl3-deficient retinal phenotype. Collectively, our study reveals an epitranscriptomic mechanism governing progenitor-to-glial cell transition during late retinogenesis, which is essential for the homeostasis of the mature retina. The mechanism revealed in this study might also apply to other nervous systems.