Spatiotemporal transcriptomic changes of human ovarian aging and the regulatory role of FOXP1

Abstract

Limited understanding exists regarding how aging impacts the cellular and molecular aspects of the human ovary. This study combines single-cell RNA sequencing and spatial transcriptomics to systematically characterize human ovarian aging. Spatiotemporal molecular signatures of the eight types of ovarian cells during aging are observed. An analysis of age-associated changes in gene expression reveals that DNA damage response may be a key biological pathway in oocyte aging. Three granulosa cells subtypes and five theca and stromal cells subtypes, as well as their spatiotemporal transcriptomics changes during aging, are identified. FOXP1 emerges as a regulator of ovarian aging, declining with age and inhibiting CDKN1A transcription. Silencing FOXP1 results in premature ovarian insufficiency in mice. These findings offer a comprehensive understanding of spatiotemporal variability in human ovarian aging, aiding the prioritization of potential diagnostic biomarkers and therapeutic strategies.

Fig. 1. Single-cell transcriptome profiling and spatial location of human ovarian cells.

a, The study flowchart. b, UMAP plots showing eight cell types (left) and age-dependent cell distribution (right). GC, granulosa cell; OO, oocyte; T&S, theca and stroma cells; SMCs, smooth muscle cells; EC, endothelial cell; MONO, monocytes. c, Representative GO terms of different cell types (left). Heat map showing top 50 marker genes in each cell type (right). P values were calculated by Fisher’s exact test. d, Each slide (i–viii) shows H&E staining, ST spot cell type predictions and characteristic markers, respectively of each cell type from left to right. The analysis was conducted in each ovarian tissue (n = 15).

Fig. 2. Cell type-specific spatiotemporal changes in transcriptional regulatory programs throughout ovarian aging.

a, Representative KEGG pathways of upregulated DEGs (top) and downregulated DEGs (bottom) compared between the O/Y, M/Y and O/M groups in eight cell types. P values were calculated by Fisher’s exact test. b, Gene set score analysis of cellular senescence pathways in various ovarian cell types of different groups, Y, M and O. Data were analyzed by two-sided Wilcoxon rank-sum tests. Box-and-whisker plots show the minimum, 25th percentile, median, 75th percentile and maximum. n = 3 per age group. c, ST spot illustrating cellular senescence and SASP score within the human ovary. The red arrow indicates follicles. d, Lipofuscin staining of different aged human ovaries. IOD, integrated optical density. Data are presented as the mean ± s.e.m. n = 10 for each group (one-way ANOVA). e, Violin plot showing the expression of CDKN1A in all ovarian cell types. Data were analyzed by two-sided Wilcoxon rank-sum tests. f, ST spot indicating the CDKN1A expression in human ovaries of different ages. g, Multiplex IHC staining for CDKN1A, ZP3 (oocyte marker), GSTA1 (GC marker) and DCN (T&S cell marker) of differently aged human ovaries. The experiment was repeated three times. h, Cytokine oligonucleotide array for SASPs in human ovaries of different ages. n = 3 for each group. Source data

Fig. 3. Spatiotemporal changes of oocytes during aging.

a, Pseudotime trajectory plot of oocytes. b, UMAP plot of oocytes from young, middle-aged and older age groups. c, Two-dimensional (2D) graph of the pseudotime-ordered oocyte cells from young, middle-aged and older age groups. d, Heat map showing the dynamic DEGs along the pseudotime. The related biological process of each subtype is listed on the right. e, Gene set score analysis of DNA damage and DNA repair pathways in oocytes of different groups. Data were analyzed by two-sided Wilcoxon rank-sum tests. Box-and-whisker plots show the minimum, 25th percentile, median, 75th percentile and maximum. n = 3 per age group. f, Cell communication analysis of oocytes between other ovarian cells. g, MDK signaling pathway network of oocytes between other ovarian cells in Y, M and O groups. h, ST spot showing the score of MDK-LRP1 in oocytes and their surround cells in ovaries of three groups, Y, M and O. The analysis was conducted in each ovarian tissue (n = 15). i, Representative image of oocytes from Y, M and O human ovary sections stained for MDK protein by IHC. The red arrows indicate primordial oocytes. The scores are listed on the right. Data are presented as the mean ± s.e.m. n = 10 for each group (one-way ANOVA). Source data

Fig. 4. Spatiotemporal transcriptional changes of three GCs subpopulations during aging.

a, UMAP visualization of GC subclusters. b, Histogram showing the cell rate of three GC subclusters in Y, M and O groups. c, Dot plot heat map showing top eight markers for subcluster GCs. d, The spatial cluster distribution of each subclusters (left) and ST spot overlapped with H&E staining (right). e, Multiplex IHC staining for the markers of subtype GCs. AMH (granulosa 1), GSTA1 (granulosa 2), LGALS1 (granulosa 3). The experiment was repeated three times. f,g, Pseudotime-ordered analysis of GCs. h, Two-dimensional graph of the pseudotime-ordered GCs from Y, M and O groups. The cell density distribution is shown in the above figure. i, Gene set score analysis of cellular senescence pathways in GC subtypes with age. Data were analyzed by two-sided Wilcoxon rank-sum tests. Box-and-whisker plots show the minimum, 25th percentile, median, 75th percentile and maximum. n = 3 per age group. j, Representative images of CDKN1A expression in GCs of three groups by IHC. k, IHC scores for CDKN1A in GCs. Data are presented as the mean ± s.e.m. (one-way ANOVA). n = 10 for each group. l, ST spot overlay of cellular senescence and SASP gene set score in GCs of three groups. m, The correlation analysis of SASP levels with age in primary human GCs (hGCs). n = 25–28 (Pearson correlation analysis, two-sided). Source data

Fig. 5. Spatiotemporal changes of five T&S cells subpopulations during aging.

a, UMAP visualization of T&S cell subclusters. b, Dot plot heat map showing top T&S cell subclusters markers. c, Heat map showing the highly expressed genes and GO terms specifically in T&S subtype cells. d, The spatial cluster distribution of each subclusters (left) and ST spot overlapped with H&E staining (right). e, IHC for markers of T&S subtype cells. The experiment was repeated three times. f, Pseudotime-ordered analysis of T&S cells (left). Subtypes are labeled by colors (right). g, A 2D graph of the pseudotime-ordered T&S cells from the young, middle and older age groups. The cell density distribution is shown in the above figure. h, Gene set score analysis of cellular senescence pathways in T&S subtype cells with age. Data were analyzed by two-sided Wilcoxon rank-sum tests. Box-and-whisker plots show the minimum, 25th percentile, median, 75th percentile and maximum. n = 3 per age group. i, Violin plot showing the expression of CDKN1A in five T&S cells subclusters. Data were analyzed by two-sided Wilcoxon rank-sum tests. j, ST spot overlay of CDKN1A expression (left) and representative images of CDKN1A expression by IHC (right) in T&S cells of three groups. k, ST spot overlay of SASPs (HMGA1, PIM1 and CAV1) expression in T&S cells of three groups. The ST data analysis is shown on the right. Data were analyzed by two-sided Wilcoxon rank-sum tests. The analysis was conducted in each ovarian tissue (n = 15). l, The correlation analysis of SASPs levels with age in pT&S cells from human ovaries. n = 25 (Pearson correlation analysis, two-sided). Source data

Fig. 6. The role of FOXP1 in ovarian aging.

a, Network of regulatory TFs in three GC subtypes. Node size correlates with the number of edges positively. b, Network of regulatory TFs in five T&S subtypes. c, The correlation analysis of FOXP1 expression with age in hGCs and pT&S. n = 28–30 (Pearson correlation analysis, two-sided). d, FOXP1 expression in ovaries. Data are presented as the mean ± s.e.m. n = 10 for each group (one-way ANOVA). e, Western blotting of FOXP1 protein levels in ovaries. Data are presented as the mean ± s.e.m. n = 6 for each group (one-way ANOVA). f, SA-β-gal staining. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). g, Protein expression of genes related to cellular senescence. The experiment was repeated for three times. h, EdU incorporation assay. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). i, Representative GO terms of DEGs. P values were calculated by Fisher’s exact test. j, The expression of SASPs and cell cycle-related genes. k, ChIP–qPCR using the FOXP1 antibody at the CDKN1A promoter. Data are presented as the mean ± s.e.m. n = 3 for each group (unpaired two-tailed t-test). l, Activity of WT and mutant (MT) CDKN1A promoter luciferase (luc) reporter. Empty vector-infected cells (MOCK) served as the control. Data are presented as the mean ± s.e.m. n = 6 for each group (one-way ANOVA). m, mRNA expression after FOXP1 overexpression in COV434. Data are presented as the mean ± s.e.m. n = 3 for each group (unpaired two-tailed t-test). n, SA-β-gal staining. The percentage of SA-β-gal-positive cells was shown on the right. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). o, mRNA expression of CDKN1A. Data are presented as the mean ± s.e.m. n = 3 for each group (unpaired two-tailed t-test). p, The protein expression of CDKN1A in COV434. Data are presented as the mean ± s.e.m. This test was repeated three times (unpaired two-tailed t-test). Source data

Fig. 7. Conditional knockout of FOXP1 in GCs accelerates ovarian aging in mice.

a, Schematic representation of the deletion of FOXP1 in GCs by using CYP19A1-Cre transgenic mice. b, Ovaries of FOXP1^+/+ and FOXP1^tr/tr mice at 12 weeks. c, H&E staining of ovaries from FOXP1^+/+ and FOXP1^tr/tr mice. d, Comparison of follicle numbers of FOXP1^+/+ and FOXP1^tr/tr mice. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). e, Serum AMH level of FOXP1^+/+ and FOXP1^tr/tr mice. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). f, Serum E2 level of FOXP1^+/+ and FOXP1^tr/tr mice. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). g, Litter size of mated mice. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). h, SA-β-gal staining of ovaries from FOXP1^+/+ and FOXP1^tr/tr mice. The number of SA-β-gal-positive cells was shown on the right. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). i, Relative RNA expression of SASPs in GCs from FOXP1^+/+ and FOXP1^tr/tr mice. Data are presented as the mean ± s.e.m. n = 3 for each group (unpaired two-tailed t-test). j, Representative images of γH2AX staining in ovaries of FOXP1^+/+ and FOXP1^tr/tr mice. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). k, TUNEL staining of ovaries from FOXP1^+/+ and FOXP1^tr/tr mice. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). Source data

Fig. 8. Quercetin treatment protects the ovarian reserve in middle-aged mice.

a, SA-β-gal staining in COV434 upon administration of fisetin (F), quercetin (Q) and dasatinib (D) in cells with knockdown of FOXP1. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). b, Relative RNA expression of CDKN1A in si-FOXP1 COV434 treated with F, Q and D. Data are presented as the mean ± s.e.m. n = 3 for each group (unpaired two-tailed t-test). c, EdU incorporation assay of COV434 upon administration of F, Q and D in cells with knockdown of FOXP1. Data are presented as the mean ± s.e.m. n = 5 for each group (unpaired two-tailed t-test). d, FOXP1 protein expression in COV434 treated by F, Q and D. This test was repeated three times. Data are presented as the mean ± s.e.m. (one-way ANOVA). n = 3 for each group. e, Experimental design to test the effects of quercetin on ovarian aging. f, Estrous cycles monitoring of mice. n = 25 for each group (Fisher Freeman Halton, two-sided). g, Serum AMH levels. Data are presented as the mean ± s.e.m. n = 8 for each group (unpaired two-tailed t-test). h, Follicle counting results according to ovary serial sections. ANF, antral follicle; ATF, atretic follicle; PMF, primordial follicle; PF, primary follicle; SF, secondary follicle; CL, corpus luteum. Data are presented as the mean ± s.d. n = 6 for each group (unpaired two-tailed t-test). i, The proportion of successful pregnant mice. n = 8 for each group (chi-squared test, two-sided). j, SA-β-gal staining of mice ovaries. k, Statistical analysis of SA-β-gal-positive GCs in mice ovaries. Data are presented as the mean ± s.e.m. (one-way ANOVA). n = 5 for each group. l, Western blot analysis of FOXP1, CDKN1A and γH2AX levels from the ovaries of treated mice. m, Densitometry quantified data of western blot analysis. Data are presented as the mean ± s.e.m. (unpaired two-tailed t-test). This test was repeated three times. n, Schematic illustration showing that downregulation of FOXP1 in GCs and T&S cells contributes to ovarian senescence. Source data

Extended Data Fig. 1. Information of human ovaries and quality control of scRNA-seq data.

(a) H&E staining of ovaries from women in young, middle and old groups for scRNA-seq and ST-seq (n = 15). (b) Percentages of mitochondrial genes detected in each sample. Box-and-whisker plots (minimum, 25th percentile, median, 75th percentile, maximum). (c) Left: numbers of unique molecular identifiers (UMIs) (upper) and genes (lower) detected in each cell. Right: UMIs counts and expressed genes are shown by three groups. (d) The number of expressed UMIs (nUMI) and genes (nGene) in different cell types. (e) UMAP plot showing the eight ovarian cell clusters in each sample. (f) UMAP plot showing all ovarian cell clusters. (g) Top three markers of each cell type shown as expression quantity (circle size) and mean expression (color). (h) UMAP cluster map of each slide showing expression of marker genes. (i) Bar plot showing the cell numbers of different cell types and detected genes in human ovaries. (j) Immunofluorescence staining of human GCs markers (AMH and GSTA1) and oocytes markers (ZP3 and TUBB8). The experiment was repeated for three times.

Extended Data Fig. 2. Quality control and summary of ST-seq data.

(a) Violin plots showing the percentage of mitochondrial genes detected in each sample. (b) Spatial feature plots of the number of expressed UMI (nUMIs) and genes (nGene) in each sample. (c) Clustered ST spots integrated with scRNA-Seq cell type annotations mapped to the H&E image showing eight cell types in each sample.

Extended Data Fig. 3. Supplement of gene expression in different cell types changed throughout human ovarian aging.

(a) Heatmaps showing the upregulated (left) and downregulated (right) DEGs of each cell type between the old and young groups (O/Y), middle and young groups (M/Y), and old and middle groups (O/M). Gene numbers on the left represents from top to bottom, DEGs shared by at least two cell types, DEGs shared by at least two groups, unique DEGs of each cell type in each group. The numbers of unique DEGs are annotated. (b) Venn diagram of DEGs shared by three groups. Y, young; M, middle; O, old. (c) Fluorescence-based-β-Gal staining of human ovaries from young, middle and old group. The scores are listed on the right. Data are presented as the mean ± SEM. n = 10 for each group (one-way ANOVA). (d) Relative mRNA expression of SASPs in human ovaries at different ages by RT–qPCR. Data are presented as the mean ± SEM. n = 9 for each group (one-way ANOVA). (e) Gene set score analysis of NF-κB signaling pathway in eight ovarian cell types of different groups. Two-sided Wilcoxon rank-sum tests. Box-and-whisker plots (minimum, 25th percentile, median, 75th percentile, maximum). n = 3 per age group. (f) Protein expression of cellular senescence-related genes in human ovaries detected by Western blot. Representative images were shown. Data are presented as the mean ± SEM. This test was repeated three times (one-way ANOVA). Source data

Extended Data Fig. 4. Supplement of gene expression in oocytes changed throughout human ovarian aging.

(a) The expression level of stage-specific markers for oocytes. LMOD3 and FOS for primordial oocytes, RPS4X and FIGLA for primary oocytes, SYT5, STK26 and TAF1A for secondary oocytes, UBOX5 and CCDC25 for antral oocytes, HTRA3 and NBPF12 for preovulotary oocytes. (b) Heat map illustrating the DEGs of each oocyte stage. DEG numbers are shown on the left. (c) Representative genes of DNA damage (STAT3, EIF4A1) and DNA repair (APEX1, RAD1) expression in oocytes of three groups. Y, young; M, middle; O, old. (d) Representative images of oocytes by IHC of γH2AX, 8-OHdG and nitrotyrosine between three groups.(e) IHC scores of relative expression for each group. AOD, average optical density. Data are presented as the mean ± SEM. n = 10 for each group (one-way ANOVA).

Extended Data Fig. 5. Supplement of cell communications in oocytes.

(a) The receptor-ligands of communication between oocytes and other ovarian cells in the old group. P values were calculated by Fisher’s exact test. (b) ST spot showing the expression of receptor–ligand pair MDK and LRP1 in human ovary.

Extended Data Fig. 6. Supplement of GCs subclusters.

(a) Heat map showing the highly expressed genes specifically in granulosa 1, granulosa 2, granulosa 3. (b) ST spot overlay of PCNA expression in GCs subtypes. (c) Heat map showing the dynamic DEGs along the pseudotime. Subtypes are labeled by colors (upper panel). The related biological process of each subtype is listed on the right. The related genes are listed on the left. (d) Heatmaps showing the upregulated (above, orange) and downregulated (below, blue) DEGs for each GCs subclusters between three groups (O/Y, M/Y, O/M). (e) Representative GO terms of DEGs enrichment between three groups (O/Y, M/Y, O/M) in each GCs subtypes. Left, the GO terms of upregulated DEGs. Right, the GO terms of downregulated DEGs. P values were calculated by Fisher’s exact test. (f) Violin plot showing the expression of CDKN1A in three GCs subclusters during aging. Two-sided Wilcoxon rank-sum tests. (g) ST spot overlay of LMNA expression in GCs of three groups. The ST data analysis is shown on the right. Two-sided Wilcoxon rank-sum tests. Box-and-whisker plots (minimum, 25th percentile, median, 75th percentile, maximum). n = 5 per age group.

Extended Data Fig. 7. Supplement of T&S cells subclusters.

(a) Histogram showing the cell rate of five T&S subclusters in young (Y), middle (M) and old (O) groups. (b) Heat map showing the dynamic DEGs along the pseudotime. Subtypes are labeled by colors (upper panel). The related biological process of each subtype is listed on the right. The related genes are listed on the left. (c) Heatmaps showing the upregulated (left, orange) and downregulated (right, blue) DEGs for each T&S subclusters between three groups (O/Y, M/Y, O/M). (d) Representative GO terms of DEGs enrichment between three groups (O/Y, M/Y, O/M) in each theca & stroma cell subtype. Left, the GO terms of upregulated DEGs. Right, the GO terms of downregulated DEGs. P values were calculated by Fisher’s exact test. (e) ST data analysis for NF-κB in T&S cells of three groups. Y, young; M, middle; O, old. Two-sided Wilcoxon rank-sum tests.

Extended Data Fig. 8. Supplement of transcription factor.

(a) Relative mRNA expression of FOXP1, SOX4, FOS-knockdown COV434 cells. Data are presented as the mean ± SEM. n = 3 for each group (unpaired two-tailed t-test). (b) Analysis of protein expression of genes related to cellular senescence in COV434 and pT&S upon si-FOXP1-mediated gene knockdown. Data are presented as the mean ± SEM. n = 3 for each group (unpaired two-tailed t-test). (c) Immunofluorescence staining of Ki67 in COV434 and pT&S upon si-FOXP1-mediated gene knockdown. The relative intensity is listed on the right. Data are presented as the mean ± SEM. n = 5 for each group (unpaired two-tailed t-test). (d) Immunofluorescence staining of γH2AX in COV434 and pT&S upon si-FOXP1-mediated gene knockdown. The relative intensity is listed on the right. Data are presented as the mean ± SEM. n = 5 for each group (unpaired two-tailed t-test). (e) PCA analysis of three si-NC samples and three si-FOXP1 samples. (f) Volcano plot of DEGs. Blue: downregulated DEGs; Red: upregulated DEGs. (g) KEGG terms associated with cellular senescence in si-NC group and si-FOXP1 group. P values were calculated by Fisher’s exact test. (h) The predicted binding sites for FOXP1 in the promoter of CDKN1A using the JASPAR website.

Extended Data Fig. 9. The effect of quercetin in COV434 with knockdown of FOXP1.

(a) mRNA relative expression of CDKN2A, IL-6 and IL-8 in si-FOXP1 COV434 treated with fisetin (F), quercetin (Q) and dasatinib (D). Data are presented as the mean ± SEM. n = 3 for each group (unpaired two-tailed t-test). (b) Immunofluorescence staining of Ki67 upon administration of fisetin (F), quercetin (Q) and dasatinib (D) in COV434 with knockdown of FOXP1. The relative intensity is listed on the right. Data are presented as the mean ± SEM. n = 5 for each group (unpaired two-tailed t-test). (c) Immunofluorescence staining of γH2AX upon administration of fisetin (F), quercetin (Q) and dasatinib (D) in COV434 with knockdown of FOXP1. The relative intensity is listed on the right. Data are presented as the mean ± SEM. n = 5 for each group (unpaired two-tailed t-test).

Extended Data Fig. 10. Biological safety evaluation of quercetin treatment in mice.

(a) The body weight of mice in four groups. Data are presented as the mean ± SEM. n = 25 for each group (unpaired two-tailed t-test). (b) and (c) The serum level of alanine aminotransferase (ALT) and aspartate transaminase (AST) for liver function. Data are presented as the mean ± SEM. n = 5 for each group (unpaired two-tailed t-test). (d) and (e) The serum level creatinine (CREA) and blood urea nitrogen (BUN) for renal function. Data are presented as the mean ± SEM. n = 5 for each group (unpaired two-tailed t-test). (f) and (g) The serum level of creatine kinase (CK) and lactic dehydrogenase (LDH) for heart function. Data are presented as the mean ± SEM. n = 5 for each group (unpaired two-tailed t-test).