Single-Cell Transcriptome Profiling Identified Eight Cell Types in the Human Ovary
We obtained ovaries from women undergoing surgeries such as hysterectomy and oophorectomy because of cervical cancer or endometrial cancer. We chose three representative groups to better understand the mechanism of ovarian aging: young group (18y, 22y, 28y), middle group (37y, 38y, 39y), and old group (47y, 48y, 49y) (Table S1). The morphological analysis revealed more primordial follicles in the young group than in the middle and old groups, but more atretic follicles were observed in the older group (Figure S1A).
We applied 10x Genomics to the ovarian single-cell analysis to further investigate the mechanisms of human ovarian aging, especially the specific characteristics of each single cell. The 0.5–1cm3 ovary samples from women undergoing surgery were enzymatically dissociated into single cells within two hours (Fig. 1A), and the cell viability was greater than 80%. We obtained 92965 ovarian cells (31005 cells from young group, 32557 cells from middle group and 29403 cells from old group). Cells expressing high levels of mitochondrial genes (> 10% of total UMIs) were excluded (Figure S1B), and after quality control, only 63612 cells were retained for further analysis. The total cellular RNA content and number of expressed genes did not differ among the three groups (Figure S1C). Then, we used the UMAP algorithm for the nonlinear dimensionality reduction analysis and identified eight cell types based on specific cell markers (Figs. 1B and S1D-S1F; Table S2). They were granulosa cells (GCs, GSTA1+, AMH+ and HSD17B1+), oocytes (OO, TUBB8+, ZP3+ and FIGLA+), theca & stroma cells (T&S, DCN+ and STAR+), smooth muscle cells (SMCs, ACTA2+ and MUSTN1+), endothelial cells (ECs, TM4SF1+ and VWF+), monocyte (MONO, TYROBP+ and IFI30+), Natural Killer cells (NK, CCL5+ and NKG7+), and T lymphocytes (T, IL7R+ and KLRB1+) (Figs. 1C and S1G). Similar to a previous study 7, 8, T&S cells accounted for the majority of ovarian cells (Figure S1H). By comparing the cell types and the percentage of each type within the three groups, we found that the cell types and most cell densities did not differ with age. However, the percentage of smooth muscle cells increased and percentage of endothelial cells decreased (no significant difference) with age (Fig. 1D), which may be associated with vascular remodeling in the microenvironment of ovarian aging. Then we validated the specificity of markers of oocytes and granulosa cells (Figure S1I). In addition to classical markers such as ZP3 for oocytes, we found that TUBB8, which encodes the primary beta-tubulin subunit expressed in oocytes 12, was also specifically distributed in oocytes. For granulosa cells, we found that GSTA1 was specifically located in GCs, indicating that GSTA1 may be a potential marker for granulosa cells. Furthermore, we performed a GO analysis of the marker genes of each cell type ensure the veracity of cell clustering (Fig. 1E). For example, the genes with high expression in GCs were enriched in the regulation of hormone levels, and those expressed at high levels in oocytes were enriched in oocyte differentiation. GO terms specific to T&S cells included ‘‘extracellular matrix organization” and “cholesterol transfer activity’’. GO terms including “muscle system process” and “cytoskeleton organization” were enriched in SMCs. Collectively, our data represent the first single-cell transcriptomic map for human ovarian aging, which illuminated the cell composition of the human ovary and the changes with age.
Spatial Location of Human Ovarian Cells
We performed ST on tissues collected from across ovarian aging (3 sections from 3 samples; young, n = 1; middle, n = 1; old, n = 1) to map spatial distributions of the cells identified using scRNA-seq data. In the human ovarian ST-seq data, the mtRNA of all ST spots represented less than 20% of the total reads, indicating good data quality (Figure S2A). Then, SCTransform was used to normalize the data (Figure S2B). ST-seq of the human ovary collectively detected over 21,000 genes (Figures S2C and S2D), and the location of follicles was detected with higher UMIs and more genes (Figure S2C).
Using the scRNA-seq atlas, we performed a factor analysis to determine the likely single-cell composition of each spot, thus spatially localizing all scRNA-seq clusters. These well-characterized cell types were localized, such as oocytes in the middle of the follicle (Fig. 2Aⅰ), granulosa cells in the outer follicle (Fig. 2Aⅱ), theca & stroma cells (Fig. 2Aⅲ) that were widely distributed, and smooth muscle cells/endothelial cells that were distributed toward the blood vessel (Figs. 2Aⅳ and 2Aⅴ). The immune cells, including monocytes, Natural Killer cells and T lymphocytes, were mainly distributed in the interstitium of the medulla, atretic follicles and corpus luteum (Figs. 2Aⅵ, 2Aⅶ and 2Aⅷ). The overall distribution of all cells is shown in Figure S2E. Furthermore, after modeling the spatially resolved expression of marker genes for different types of cells, we observed that the local expression of the oocyte marker genes ZP3 and TUBB8 predicted the location of oocytes (Fig. 2Aⅰ), the granulosa cells marker genes AMH and GSTA1 predicted the location of granulosa cells (Fig. 2Aⅱ), the theca & stroma cells marker genes DCN and STAR predicted the location of theca & stroma cells (Fig. 2Aⅲ), ACTA2 and MUSTN1 predicted the location of smooth muscle cells (Fig. 2Aⅳ), TM4SF1 and VWF predicted the location of endothelial cells (Fig. 2Aⅴ), TYROBP and IFI30 predicted the location of monocyte (Fig. 2Aⅵ), NKG7 predicted the location of Natural Killer cells (Fig. 2Aⅶ), IL7R predicted the location of T lymphocytes (Fig. 2Aⅷ).
The human ovarian structure includes the cortex, transition region and medulla. We wanted to determine whether the transcriptional/cellular spatial variability corresponds to the tissue depth. We identified 142 depth-associated genes, reflecting pathways active in different layers: (Fig. 2B)—the surface region enriched for the response to temperature stimulus, regulation of response to stimulus and response to incorrect protein progressing toward response and stimulation to the outside, indicating that the ovarian surface is vulnerable to internal and external stimulation; and, deeper spots enriched for translational initiation, RNA catabolic process, and cotranslational protein targeting to membrane progressing toward gene expression, indicating that there is an active biological process in the deeper part of the ovaries.
Gene Expression in Different Cell Types Changed Throughout Human Ovarian Aging
We compared the gene expression patterns of ovarian cell types between the young, middle and old groups to further explore the mechanism of ovarian aging at the cellular level. We identified thousands of differentially expressed genes (DEGs, |avg_logFC| > 0.25 and p_val_adj < 0.05) in at least one cell type of the human ovary during aging (Figs. 3A and 3B; Table S4). A total of 1068, 711, and 889 upregulated DEGs were identified between the young and old (O/Y) groups, young and middle-aged (M/Y) groups, and middle-aged and old (O/M) groups, respectively (Fig. 3C). In comparison, 1187, 376, and 1241 downregulated DEGs were identified in the O/Y groups, M/Y groups, and O/M groups, respectively (Fig. 3C). Among all ovarian cell types, oocytes showed the greatest difference among the three groups (Figure S3A). Notably, the analysis of DEGs revealed a substantial difference between the perimenopausal ovaries and the young or middle-aged ovaries with reproductive function.
We aligned datasets from individual samples by their chronological age and clustered the DEGs by their expression patterns to identify the constantly upregulated or downregulated DEGs during aging. Using this method, we identified 5,476 age-dependent upregulated DEGs and 3,412 age-dependent downregulated DEGs in eight cell types (Fig. 3D; Table S5). GO and KEGG enrichment analyses revealed that age-related upregulated genes were mainly associated with cellular senescence, the FoxO signaling pathway, the IL-17 signaling pathway, the nuclear factor-κB (NF-κB) signaling pathway, the NOD-like receptor signaling pathway, the p53 signaling pathway, the PI3K-Akt signaling pathway and the transcriptional misregulation in cancer (Fig. 3E). Age-related downregulated genes were mainly related to cell migration, ECM-receptor interaction, estrogen signaling pathway, extracellular vesicle, oxidative phosphorylation, platelet activation, regulation of actin cytoskeleton and tight junction (Fig. 3E). Similarly, a pathway analysis of the spatial gene expression data indicated increased enrichment of cellular senescence and the p53 signaling pathway, while we identified decreased enrichment of estrogen signaling pathway and oxidation phosphorylation within the old ovary (Fig. 3F).
Based on the aforementioned results, cellular senescence, a process that results from a variety of stresses and leads to a state of irreversible growth arrest, may play an important role in human ovarian aging. We compared the core senescence gene list from the GeneAge database to our scRNA-seq data and generated a cellular senescence score to assess the senescent cell fraction 13. The cellular senescence score was significantly increased in multiple cell types during ovarian aging (Figs. 3G and S3B). We further detected lipofuscin, which accumulates mainly in aged cells and is considered a hallmark of cellular senescence 14, in human ovaries from patients of different ages. As shown in Fig. 3H, lipofuscin accumulation was increased in aged ovarian tissues. The expression of CDKN1A/p21, the core senescence marker, was significantly increased in ovarian cells from the old group compared to the young or middle group (Figure S3C). Additionally, we identified the spatial expression of CDKN1A in the three groups using spatial transcriptomic data, and showed increased activity in the ovarian cortical region during aging (Figure S3D). We further validated the finding by performing immunostaining and found that increased numbers of CDKN1A-positive cells accumulated during ovarian aging (Fig. 3I). Consistently, we also observed increased CDKN1A mRNA expression in aged ovaries (Figure S3E), indicating that increased numbers of senescent cells may underlie ovarian aging. Furthermore, cellular senescence is also accompanied by chronic inflammation. We observed that inflammatory response genes and the NF-kB signaling pathway, which in turn induce the senescence-associated secretory phenotype (SASP) and aggravate cellular senescence 15, were activated in most cell types during ovarian aging (Figures S3F and S3G). Collectively, these results indicate that a series of molecular changes, especially cellular senescence and the inflammatory response, occur in aged human ovaries and likely contribute to the development of aging-related ovarian disorders.
Changes in Oocytes during Human Ovarian Aging
We next sought to identify aging-associated changes in gene expression in oocytes. We performed a pseudotime trajectory analysis to analyze the origin and maturation of oocytes (Fig. 4A). All oocytes in the young group were identified as Stage_a and Stage_b cells, most oocytes in the middle group were identified as Stage_a and all oocytes in the old group were predominantly classified in Stage_c cells (Fig. 4B). This pseudotime analysis showed that the oocyte in young group were at the beginning of the trajectory path, whereas the oocytes in old group were at a terminal state, which showed that the oocytes undergo GeneSwitch in the process of aging (Fig. 4C). Next, we further analyzed the dynamic changes in the gene expression pattern at each stage. A total of 356, 219, and 3 DEGs (p adj_val < 0.05, log2FC > 0.5) were observed in the three states, respectively (Fig. 4D). Then, we performed a GO analysis of DEGs at each stage. In Stage_a cells, the main enriched GO terms were “cotranslational protein targeting to membrane,” “nuclear-transcribed mRNA catabolic process,” and “translational initiation” (Fig. 4E), suggesting that oocytes at this stage exhibit active biological processes. In Stage_b cells, “collagen-containing extracellular matrix”, “extracellular matrix organization”, and “positive regulation of cell adhesion” were enriched (Fig. 4E). In Stage_c cells, the main enriched GO terms were “organelle fission”, “polysomal ribosome” and “double-strand break repair” (Fig. 4E), indicating that oocytes in Stage_c undergo DNA damage and cell stress (Fig. 4E). A recent study revealed that the DNA damage response (DDR) is the primary biological pathway that regulates reproductive senescence 16. We compared the DNA damage and repair gene list and generated DNA damage and repair scores. With aging, DNA damage genes were upregulated in oocytes (Fig. 4F). The expression of classical DNA damage response genes, such as STAT3 and EIF4A1, was markedly upregulated in old oocytes (Fig. 4F). However, DNA repair genes that are crucial for the maintenance of oocyte homeostasis, such as APEX1 and RAD1, were downregulated in old oocytes (Fig. 4G). Consistently, we observed the accumulation of DNA oxidation (8-OHdG-positive cells), DNA damage (γH2AX-positive cells), and protein oxidation (nitrotyrosine-positive cells) in aged oocytes (Figs. 4H and 4I).
Changes in the Transcriptional Profiles of the Three Subpopulations of Granulosa Cells during Human Ovarian Aging
Granulosa cells of the ovarian follicle surround and interact with the developing oocyte, and are responsible for estrogen and progesterone synthesis. First, we explored the heterogeneity of human granulosa cells. UMAP analysis classified all the granulosa cells into three groups (granulosa cell subtypes 1–3) (Fig. 5A; Table S2). Statistics indicate that granulosa cells in the young group were mainly granulosa cell subtype 1 and 2, and the granulosa cells in the old group were mainly granulosa cell subtype 3 (Fig. 5B). Further analysis of the expression levels of the markers showed granulosa cell subtype 1 was characterized by the expression of the known markers AMH, FST, HSD17B1, SERPINE2, and PRKAR2B, as well as several genes that had not previously been associated with granulosa cells, such as TNNI3, DSP, and MAGED2 (Fig. 5C). The GO analysis revealed the biological processes that were enriched for follicular development. Notably, ‘‘ATP metabolic process’’ and “Gap junction’’ were enriched in granulosa cell subtype 1 (Fig. 5D). Granulosa cell subtype 2 was characterized by the expression of genes related to hormone synthesis, such as INSL3, APOE, GSTA1, APOA1, FDX1 and CYP17A1 (Fig. 5C). Similarly, GO terms, including ‘‘cholesterol metabolic process’’ and ‘‘steroid biosynthetic process’’ were enriched in granulosa cell subtype 2 (Fig. 5D). In addition, granulosa cell subtype 3 by the expression of the markers DCN, LGALS1 and LGALS1, which were reported to regulate the apoptosis and cell cycle of granulosa cells 17, 18 (Fig. 5C). Notably, ‘‘Apoptosis’’ and ‘‘Cell cycle’’ were enriched in granulosa cell subtype 3 (Fig. 5D).
Since each individual spot on the spatial transcriptomic slide putatively contains multiple cells, we transferred the labels from the integrated data to spatial gene expression data and mapped different cell types based on spatial location (Fig. 5E). Interestingly, we observed three different populations of granulosa cells in three distinct areas. Granulosa cell subtype 1 was located in the cumulus of the antral follicle, granulosa cell subtype 2 was located in mural layer of follicles, while granulosa cell subtype 3 showed a broad distribution across antral follicles and the ovarian cortex (Fig. 5E).
We next explored the dynamic states and cell transitions in granulosa cells by inferring the state trajectories using Monocle. This analysis showed that granulosa cell subtype 1 was at the beginning of the trajectory path (Figs. 5F and 5G). We analyzed the trajectories of granulosa cells in the three groups separately to further delineate the transition states associated with granulosa cells in different samples. Surprisingly, early-stage granulosa cell subtype 1 was predominantly distributed in young samples, whereas granulosa cells in old samples were primarily at the terminal ends of the granulosa cell subtype 3, indicating that the granulosa cells undergo GeneSwitch during the ovarian aging (Fig. 5H). We next investigated the transcriptional changes associated with transitional states and observed that the granulosa cell clusters were categorized into 3 phases (Fig. 5I). Granulosa cell subtype 1 was predominantly phase 1 cells, and pathway analysis indicated that signaling pathways involved in the metabolism of GTP and glucose. This finding indicated that phase 1 cells were more closely related to metabolism and that these metabolites may reach the oocyte via paracrine signals and gap junctions to promote follicular development. Pathway analysis suggested that cells in phase 2 were enriched in the apoptosis pathway, matching the characteristics of granulosa cell subtype 3. Phase 3 was characterized by genes involved in cytokine production and extracellular matrix, further confirming the hormone synthesis function of granulosa cell subtype 2. Collectively, the above results indicated that granulosa cells can be divided into three subtypes, with granulosa cell subtypes 1 and 2 have better function, mainly in the young and middle ovaries, while granulosa cell subtype 3 was mainly present in the old ovaries.
As apoptosis of granulosa cells causes follicular atresia and ovarian aging (Matsuda et al., 2012), we next focused on the aging-associated changes in gene expression in granulosa cells. Further dissection of the transcriptomic changes in these subpopulations during aging showed that DEGs in granulosa cells between the young, middle and old groups partially overlapped among different subpopulations (Figures S4A and S4B). We also identified age-related DEGs in each subpopulation by performing a hierarchical clustering analysis (Figure S4C; Table S5). We calculated the cellular senescence score to determine whether cell senescence are also changed in granulosa cells during aging, and the results showed that it was augmented in all types of granulosa cells in the aged ovary (Fig. 5J). Accordingly, CDKN1A, a molecular marker of senescent cells, was upregulated in three subpopulations of granulosa cells (Figs. 5K and S4D). We obtained human granulosa cells (hGCs) from healthy women aged 21 to 46 years old to further validate the senescence of granulosa cells (Figures S4E). We also observed aging-associated increases in the expression CDKN1A, CDKN2A and SASP (IL-6 and IL-8) in hGCs from healthy donors, further supporting cellular senescence in hGCs during physiological aging. Finally, we performed a transcriptional regulatory network analysis for each granulosa cell subpopulation and identified preferential transcriptional regulons in each subpopulation (Fig. 5L). FOS, SOX4, FOXP1 and KLF2 served as core downregulated transcription factor (TF), while JUND, CEBPD, FOSB and IRF1 served as core upregulated transcription factor at the hub of transcriptional regulatory networks of granulosa cells in different subpopulations during aging (Figs. 5L and S4F). Altogether, these findings suggest that cellular senescence is a feature of aging granulosa cells, likely contributing to ovarian dysfunction during aging.
Changes in the Transcriptional Profiles of the Five Subpopulations of Theca & Stroma Cells during Human Ovarian Aging
The ovary contains a large stromal compartment, including the tunica albuginea, interstitial stromal and theca interna 19. However, the classification of ovarian stromal cells and whether stromal cell subpopulations respond differently to aging are not clear. Using unbiased clustering and UMAP analyses, we identified five types of theca & stroma cells with distinct cellular transcriptomic signatures (Fig. 6A; Table S2). Furthermore, we analyzed the changes in the gene expression pattern in each subpopulation of theca & stroma cells. An analysis of the expression levels of the markers showed that theca & stroma cells subtype 1 was characterized by the expression of the markers STAR and CYB5A (Fig. 6B). The GO analysis revealed that the biological processes RNA and protein synthesis were enriched in theca & stroma cell subtype 1 (Fig. 6C). Theca & stroma cell subtype 2 was characterized by the expression of genes associated with the extracellular matrix, such as COL1A1, COL3A1 and COL1A2 (Fig. 6B). Accordingly, the GO analysis revealed that “Extracellular vesicle” and “Extracellular matrix” were enriched in theca & stroma cell subtype 2 (Fig. 6C), consistent with their function in maintaining the morphology and function of follicles. Theca & stroma cell subtype 3 expressed markers of myofibroblast-like fibroblast states, as evidenced by the increased expression of ACTA2 and TAGLN (Fig. 6B). These cells were further enriched for pathways related to supramolecular fiber, muscle system process and cytoskeletal protein binding (Fig. 6C). Theca & stroma cell subtype 4 was defined by FBLN1 and CXCL2 expression (Fig. 6B). Consistent with gene expression patterns, the GO analysis revealed an enrichment of “response to chemical” and “regulation of response to stimulus” (Fig. 6C). Theca & stroma subtype cell 5 had features of inflammatory-like fibroblasts, with high expression of markers (CD74, HLA-DRB1, CCL4 and HLA-DRA) and pathways related to “antigen processing and presentation” and “cellular response to interferon-gamma and complement cascades” (Figs. 6B and 6C). Statistics indicate decreased numbers of theca & stroma cell subtype 2 and theca & stroma cell subtype 5 in the old group, consistent with the decrease in the number of follicles in the old ovary (Figure S5A).
We spatially located five subpopulations of theca & stroma cells to integrate the scRNA-seq and ST data (Fig. 6D). Interestingly, we observed the theca & stroma cell subtype 1 was widely distributed in the ovarian medulla. Theca & stroma cell subtype 2 was distributed around the follicle, and may be involved in the development of follicles. Theca & stromal cell subtype 3 was located around the blood vessels, and the GO analysis showed that these cells are probably a group of stroma cells that maintain vasoconstriction. Theca & stroma subtype cell 4 was distributed in the outer cortex, consistent with the GO analysis, and play roles in ovarian defenses against external stress. Theca & stroma cell subtype 5 was distributed in the repair area after ovulation, consistent with their immunoregulatory properties, suggesting that they are involved in the inflammatory response and ovarian repair after ovulation. Taken together, the results highlighted that the spatiotemporal dynamic landscape of theca & stromal cells subpopulations were closely related to their function.
Furthermore, we selected the regions of different follicle stages for comparison to obtain insights into the spatial organization of theca & stroma cell types. According to H&E staining of the ovary section, region 1 represents healthy follicle, and the regions 2 and 3 represent atretic follicles. Region 4 is the disappearing corpus luteum after ovulation with many monocytes and endothelial cells. In region 5, theca & stroma cells have replaced the corpus luteum (Figure S5B). The distribution of different types of theca & stroma cells was associated with their appropriate function (Figure S5C), and they showed distinct gene expression patterns (Figure S5D). By comparing the transcriptome datasets of healthy follicle (region 1) with that of the atretic follicle (region 2), we observed that the process of “ferroptosis, lysosome, autophagy, apoptosis and necroptosis” was enriched in region 2 (Figure S5E), suggesting that multiple cell death patterns are involved in the regulation of follicular atresia. We also analyzed the changes of the absorption process of the corpus luteum (region 4 and region 5), the data demonstrated theca & stroma cells in region 4 were enriched in the process of “ferritin complex, sequestering of iron ion, ferric iron binding, ferroptosis and mineral absorption” (Figure S5E), highlighting that iron metabolism plays an important role in the absorption of the corpus luteum.
T-SNE analysis helped to reveal heterogeneity among theca & stroma cells, and we also wondered if they shared common differentiation trajectories. Theca & stroma cells from different subclusters were distributed broadly across the pseudotime space, with theca & stroma cell subtype 2 primarily occupying the beginning of the trajectory path, whereas the remaining half mainly consisted of theca & stroma cell subtype 1 (Fig. 6E). We analyzed the trajectories of theca & stroma cells in young, middle and old groups separately to further delineate the transition states associated with theca & stromal cells in samples from patients of different ages. Surprisingly, early-stage theca & stroma cells were predominantly distributed in young and middle samples, with few cells identified at the end of the cell state transition path, whereas theca & stroma cells in old samples were primarily located at the terminal ends of the transition path (Fig. 6F). We next investigated the transcriptional changes associated with transitional states and observed that the theca & stromal cell clusters were categorized into 3 phases (Fig. 6G). Theca & stroma cell subtype 2 was predominantly phase 1 cells involved in the metabolic process and extracellular matrix, suggesting that these cells play an important role in follicle development. The other types of theca & stroma cells were all distributed in phase 2 and phase 3, suggesting limit changes in differentiation trajectories.
Furthermore, we analyzed age-dependent alterations in gene expression in the theca & stroma cells represented in our dataset. Dissection of the transcriptomic changes in these subpopulations during aging showed that DEGs in theca & stroma cells between the three groups largely overlapped among different subpopulations (Figure S5F). We also identified age-related DEGs in each subpopulation by performing a hierarchical clustering analysis (Figure S5G; Table S5), which again showed that many of them were shared across subpopulations. In addition, we observed higher cell senescence scores in all subpopulations of theca & stroma cells during aging (Fig. 6H). Accordingly, CDKN1A was upregulated in aged theca & stroma cells, which may be a consequence of ongoing cellular senescence (Figs. 6I and S5H). Finally, we performed a transcriptional regulatory network analysis for each theca & stroma cells subpopulation and identified preferential transcriptional regulons in each subpopulation (Fig. 6J). MYC, CEBPD, CREM and IRF1 were activated in a cell type-specific manner during ovarian aging, serving as core upregulated transcription factor. FOXP1, FOS, SOX4, ARX and JUN (Figure S5I), which showed decreased expression in aged theca & stroma cells, served as core downregulated transcription factor.
FOXP1 Exerted Geroprotective Effects on Human Ovarian Aging
Based on the results described above, cellular senescence may be involved in the regulation of ovarian aging. The transcription factor analysis of granulosa cells and theca & stroma cells indicated that FOXP1, SOX4 and FOS, which regulate the cellular senescence of both granulosa cells and theca & stroma cells, may be the core TFs involved in ovarian aging. We knocked down FOXP1, SOX4 and FOS using siRNAs in a human granulosa cell line (COV434) and primary theca & stroma cells (pT&S) to determine the roles of FOXP1, SOX4 and FOS in ovarian aging. Inspiringly, si-FOXP1 and si-SOX4 transfection in COV434 and pT&S cells resulted in a stable increase in the expression of cell senescence-associated genes, including CDKN1A and CDKN2A, but not si-FOS group (Figure S6A). We further examined the expression of FOXP1 and SOX4 in human ovaries and observed decreased expression of both with age (Figs. 7A and 7B). Similarly, we observed aging-associated increases in FOXP1 and SOX4 expression in hGCs from healthy donors (Fig. 7C). We speculate that the transcription factor FOXP1 and SOX4 may affect ovarian function by regulating cell senescence. Then, we silenced FOXP1 or SOX4 expression by using small interfering RNAs (siRNAs) in COV434 and pT&S cells. SA-β-gal staining showed that FOXP1 knockdown promoted the senescence of COV434 and pT&S cells, and si-SOX4 only lead to the senescence of COV434 cells (Figs. 7D and S6B). In addition, the protein expression of cellular senescence-, inflammation- and DNA damage-related genes increased upon knockdown of FOXP1 and SOX4 in COV434 or pT&S cells (Figs. 7E and S6C). Moreover, FOXP1 and SOX4 knockdown decreased the percentages of EdU-positive and Ki67-positive cells (Figs. 7F and S6D-6G), indicating that the downregulation of FOXP1 and SOX4 in human ovarian granulosa cells and theca & stroma cells led to decreased proliferation. The activation of the DNA damage response is one of the main mediators of cell senescence, which results in irreversible growth arrest 20. We also observed an increased fluorescence intensity of γH2AX in COV434 and pT&S cells upon knockdown of FOXP1 and SOX4 (Figures S6D, S6E and S6H). Together, these data indicate that an age-related decrease in FOXP1 and SOX4, especially FOXP1, which regulates the senescence of both granulosa cells and theca & stroma cells, may be an important contributor to ovarian aging.
We further explored the molecular mechanism underlying the age-related downregulation of FOXP1 and its role in human ovarian aging. Genome-wide RNA-seq analysis further revealed that FOXP1 knockdown in COV434 cells upregulated 118 genes and downregulated 405 genes (Figures S6I-S6K). The GO analysis showed that the upregulated genes were enriched in the terms ‘‘response to stimulus’’ and ‘‘immune response’’ upon knockdown of FOXP1, whereas the downregulated genes were mainly associated with “cellular developmental process” and “cell differentiation” (Fig. 7G). At the same time, the expression of senescence markers (CDKN1A, CDKN2A, and CDKN2B) the SASP (IL6, TNF, and CXCL8) was significantly higher in the si-FOXP1 groups than in the controls (Fig. 7H). Moreover, si-FOXP1 resulted in an increase in the expression of molecules in the NF-κB, TNF, NOD-like receptor, p53 and DNA damage pathways, which are both associated with cell senescence and immune inflammation (Figure S6L). We next asked whether these transcriptional changes in the human granulosa cell line were similar to those in human GCs during ovarian aging. Fifteen genes were upregulated both in FOXP1-knockdown COV434 cells and GCs from aged human ovaries, including IFITM1, IL1R1, ISG20, CFH, SOD2, CXCL2 and ISG15, which are involved in immune inflammation homeostasis (Figure S6M). In addition, ChIP-PCR was performed using the FOXP1 antibody in COV434 cells to verify whether FOXP1 directly bind to the CDKN1A promoter (Fig. 7I). These findings suggested that FOXP1 restrains ovarian cell senescence by repressing CDKN1A transcription during ovarian aging.
We subsequently explored whether three well-studied senolytics (fisetin, quercetin and dasatinib) fighting cell senescence reversed ovarian aging 21. Encouragingly, we observed that quercetin and fisetin delayed cellular senescence induced by FOXP1 knockdown (Figs. 7J and 7K). Accordingly, we observed reduced levels of senescence and inflammatory markers in fisetin- or quercetin-treated COV434 cells (Fig. 7L). In addition, fisetin and quercetin improved the proliferation rate of COV434 cells (Figs. 7M and S7A). Immunofluorescence staining also showed reduced levels of the DNA damage response marker γ-H2AX after treatment with fisetin or quercetin (Figure S7B). On the other hand, quercetin significantly decreased the levels of the CDKN1A and CDKN2A mRNAs, and both fisetin and quercetin decreased the secretion of senescence-associated cytokines (Figure S7C). These findings suggested that senolytics may prevent ovarian cells from senescence and are promising senolytics for ovarian aging, especially quercetin.
Vascular Remodeling during Human Ovarian Aging
The ovarian vasculature is important for follicle activation, growth, survival and the generation of a corpus luteum capable of maintaining pregnancy 22, while little is known about the changes in blood vessels during human ovarian aging. Morphologically, senile blood vessels exhibited characteristics of vascular aging, including increased wall thickness, increased fibrosis and lipid accumulation of the tube wall, determined by histological analysis, Masson’s staining, Sirius red staining and Oil red staining (Figs. 8A and S8A). As verified by immunofluorescence staining, we observed that CD31, a marker of ECs, labeled a curly shape in aged ovaries (Fig. 8B). Furthermore, the expression of α-SMA, a characteristic marker of SMCs, was upregulated in the vessels of aged ovaries (Fig. 8B). Together, the age-dependent structural changes in the vascular system in the human ovary include the endothelium and vascular smooth muscle cells.
We further explored age-associated ovarian vascular transcriptomics at single-cell resolution. The cells from vessels were grouped into five clusters, including endothelial cells (arterial ECs, venous ECs and capillary ECs) and smooth muscle cells (arterial SMCs and venous SMCs) (Fig. 8C). ECs were present at a higher percentage in the young group than in the middle or old group, whereas the percentages of arterial SMCs in the old group were higher than those in the young or middle group (Figure S8B). Next, scRNA-seq and ST datasets were integrated to characterize the locations of ECs and SMCs. We observed that the distribution of ECs and SMCs overlapped with the location of blood vessels (Figure S8C), further confirming the accuracy of cell clustering.
Annotations were guided by the expression of canonical cell class markers, as exemplified in Figs. 8D and S8D. Consistent with the function of these marker genes, GO and KEGG analysis revealed unique transcriptional features and enriched pathways relevant to their distinct physical functions (Figs. 8E, S8E and S8F). For example, arterial ECs were primarily associated with the KEGG term “intercellular communication”. Notably, venous ECs exhibited a predominance of vascular remodeling pathways. In addition, ‘‘extracellular matrix’’ and ‘‘blood vessel morphogenesis’’ were enriched in capillary ECs. We next investigated the gene expression patterns in SMCs, and the results showed that muscle contraction and extracellular matrix were enriched in arterial SMCs and venous SMCs, respectively (Figs. 8E and S8F).
We next compared the single-cell transcriptomes of ECs and SMCs during ovarian aging. We compared cell-intrinsic gene expression programs of the same cell types between the three groups and identified 778 upregulated and 522 downregulated genes (log2 [fold change] ≥ 0.5 and adjusted p value ≤ 0.05, referred to as young/middle/old differentially expressed genes), that were differentially expressed in at least one type of vascular cell (Fig. 8F). The core pathways annotated for upregulated DEGs that were shared by diverse vascular cell types involved inflammation and aging pathways, such as response to stress, immune response, cell death and cell senescence (Fig. 8G). In comparison, annotated pathways for downregulated DEGs that were shared by disparate vascular cell types included vasculature development, angiogenesis, extracellular matrix and estrogen signaling pathway (Fig. 8G). These observations were consistent with vessel abnormalities that we observed in aged human ovarian arteries.
According to the results shown in Fig. 3, ECs and SMCs from human ovarian tissue exhibit characteristic transcriptomic features of cellular senescence during aging. Cellular senescence may result in proatherosclerotic, proinflammatory, and prothrombotic changes in vascular function and contribute to the age-related vascular disease. Furthermore, we evaluated whether the senescence of ECs and SMCs drives vascular aging. The expression of the core senescence marker CDKN1A was significantly increased in all subtypes of ECs and SMCs (Fig. 8H). Similarly, we observed increased levels of the CDKN1A protein during ovarian aging (Fig. 8I). Moreover, the level of SASP score was significantly increased in ECs and SMCs, indicative of the cellular senescence state in the aged ovary (Fig. 8J). Furthermore, the cellular senescence associated pathway, NOD-like receptor signaling pathway and TNF signaling pathway were increased in ECs and SMCs during ovarian aging (Figure S8G). Similarly, the TF FOXP1 was decreased in the SMCs of middle and old group compared with that of young group (Figure S8H). Cellular senescence promote age-related extracellular matrix remodeling 23. SMCs are the major source of ECM, especially collagens and proteoglycans. Noticeably, the expression of Timp1, a major contributor to vascular remodel, was increased in both arterial and venous MSCs (Figure S8I). Altogether, ECM-senescent cells interactions might contribute to vascular remodel in aged human ovaries.