Identification of Key Genes and Potential New Biomarkers for Ovarian Aging: A Study Based on RNA-Sequencing Data

Abstract

Ovarian aging leads to reproductive and endocrine dysfunction, causing the disorder of multiple organs in the body and even declined quality of offspring's health. However, few studies have investigated the changes in gene expression profile in the ovarian aging process. Here, we applied integrated bioinformatics to screen, identify, and validate the critical pathogenic genes involved in ovarian aging and uncover potential molecular mechanisms. The expression profiles of GSE84078 were downloaded from the Gene Expression Omnibus (GEO) database, which included the data from ovarian samples of 10 normal C57BL/6 mice, including old (21-22 months old, ovarian failure period) and young (5-6 months old, reproductive bloom period) ovaries. First, we filtered 931 differentially expressed genes (DEGs), including 876 upregulated and 55 downregulated genes through comparison between ovarian expression data from old and young mice. Functional enrichment analysis showed that biological functions of DEGs were primarily immune response regulation, cell-cell adhesion, and phagosome pathway. The most closely related genes among DEGs (Tyrobp, Rac2, Cd14, Zap70, Lcp2, Itgb2, H2-Ab1, and Fcer1g) were identified by constructing a protein-protein interaction (PPI) network and consequently verified using mRNA and protein quantitative detection. Finally, the immune cell infiltration in the ovarian aging process was also evaluated by applying CIBERSORT, and a correlation analysis between hub genes and immune cell type was also performed. The results suggested that plasma cells and naïve CD4⁺ T cells may participate in ovarian aging. The hub genes were positively correlated with memory B cells, plasma cells, M1 macrophages, Th17 cells, and immature dendritic cells. In conclusion, this study indicates that screening for DEGs and pathways in ovarian aging using bioinformatic analysis could provide potential clues for researchers to unveil the molecular mechanism underlying ovarian aging. These results could be of clinical significance and provide effective molecular targets for the treatment of ovarian aging.

Keywords: GEO database; bioinformatics; biomarker; immune cell infiltration; ovarian aging.

FIGURE 1

Boxplot, volcano plot, and heatmap of gene expression of the dataset GSE84078. Box and whisker plot of GSE84078 before normalization (A) and after normalization (B). Volcano plot of old and young groups (C). Heatmap for selected differentially expressed genes between old and young ovary groups. A total of 931 DEGs are screened in normal old and young mice ovaries. A process of gradual color change from red to blue indicates expression values changing from high to low (D). DEGs, differentially expressed genes; FC, fold change.

FIGURE 2

Functional enrichment analysis of DEG. GO enrichment analysis of upregulated DEG in ovarian aging (A). GO enrichment analysis of downregulated DEG in ovarian aging. GO analysis divided DEG into three functional groups: biological process (BP), cellular component (CC), and molecular function (MF) (B). KEGG pathway enrichment analysis of DEG (C). Distribution of DEG in ovarian aging in the significant GO terms of biological process (D). Chord plot of the relationship between DEG and KEGG pathways (E).

FIGURE 3

PPI network complex of DEGs in ovarian aging. Red nodes represent upregulated genes; green nodes represent downregulated genes. Round rectangle nodes represent hub genes (A). Four of the highly connected clusters were identified by centrality calculating algorithms and visualized by the Centiscape plugin of Cytoscape. Interactions are color coded according to combined scores with darker edges corresponding to higher scores (B). Plot of hub genes displayed in top GO terms and KEGG pathways. Size of circles represents the relative mRNA expression of the genes (C).

FIGURE 4

Gene set enrichment analysis of GO terms (A). Gene set enrichment analysis of KEGG pathways (B). Gene set variation analysis and the differentially expressed analysis of pathways (C).

FIGURE 5

Immune cell infiltration evaluation and correlation analysis. The composition of 25 immune cell types in each sample (A). Correlation heatmap of the immune cell types (B); the intensity of the color indicated the strength of the correlation; red represents a positive correlation, while blue represents a negative correlation; violin plot of the proportion of various immune cell types (C); p value < 0.05 was marked red. Correlation between Lcp2 and infiltrating immune cells (D). Correlation between H2-Ab1 and infiltrating immune cells (E).

FIGURE 6

Correlation analysis of Lcp2, H2-Ab1, and their responding immune cell types. Correlation analysis of Lcp2 and the listed immune cell types (A): memory B cells (a), plasma cells (b), m1 macrophages (c), Th17 cells (d), and immature dendritic cells (e). Correlation analysis of H2-Ab1 and the listed immune cell types (B): memory B cells (a), plasma cells (b), m1 macrophages (c), Th17 cells (d), immature dendritic cells (e), naïve CD4⁺ T cells (f), Th2 cells (g), and resting NK cells (h).

FIGURE 7

Experiment validation of hub genes with top high degree scores. Relative mRNA expression levels of hub genes expression in the 24-week-old and 1.5-year-old mice ovary (unpaired Student’s t test) (A). Representative images of IHC detection of selected hub genes (Lcp2, H2-Ab1, Itgb2, and CD14) expression in the 24-week-old and 1.5-year-old mice ovary (B). Relative expression level analysis of LCP2, IGTB2, and CD14 (unpaired Student’s t test) (C). *p < 0.05; **p < 0.005; ***p < 0.001; ****p < 0.0001.