Abstract Background Premature ovarian insufficiency (POI) is a reproductive disorder characterized by the cessation of ovarian function before the age of 40. Although mitochondrial dysfunction and immune disorders are believed to contribute to ovarian damage in POI, the interplay between these factors remains understudied. Methods In this research, transcriptomic data related to POI were obtained from the NCBI GEO database. Hub biomarkers were identified through the construction of a protein‒protein interaction (PPI) network and further validated using RT‒qPCR and Western blot. Moreover, their expression across various cell types was elucidated via single-cell RNA sequencing analysis. A comprehensive investigation of the mitochondrial and immune profiles of POI was carried out through correlation analysis. Furthermore, potential therapeutic agents were predicted utilizing the cMap database. Results A total of 119 mitochondria-related differentially expressed genes (MitoDEGs) were identified and shown to be significantly enriched in metabolic pathways. Among these genes, Hadhb, Cpt1a, Mrpl12, and Mrps7 were confirmed both in a POI model and in human granulosa cells (GCs), where they were found to accumulate in GCs and theca cells. Immune analysis revealed variations in macrophages, monocytes, and 15 other immune cell types between the POI and control groups. Notably, strong correlations were observed between seven hub-MitoDEGs (Hadhb, Cpt1a, Cpt2, Mrpl12, Mrps7, Mrpl51, and Eci1) and various functions, such as mitochondrial respiratory complexes, dynamics, mitophagy, mitochondrial metabolism, immune-related genes, and immunocytes. Additionally, nine potential drugs (calyculin, amodiaquine, eudesmic acid, cefotaxime, BX-912, prostratin, SCH-79797, HU-211, and pizotifen) targeting key genes were identified. Conclusions Our results highlight the crosstalk between mitochondrial function and the immune response in the development of POI. The identification of MitoDEGs could lead to reliable biomarkers for the early diagnosis, monitoring, and personalized treatment of POI. Supplementary Information The online version contains supplementary material available at 10.1186/s12916-024-03675-7. Keywords: Premature ovarian insufficiency, Mitochondrial dysfunction, Immune cell infiltration, Bioinformatics analysis Background Premature ovarian insufficiency (POI) is a significant condition characterized by the depletion of ovarian function before the age of 40, affecting at least 1% of women of reproductive age [[36]1]. Individuals with POI often experience subfertility and symptoms commonly associated with menopause, such as hot flashes, night sweats, and vaginal dryness, which can greatly impact their quality of life [[37]2, [38]3]. Unfortunately, existing treatments are limited in effectiveness due to the incomplete restoration of ovarian function and the presence of severe side effects [[39]4, [40]5]. Consequently, there is an urgent need for advancements in the diagnosis and management of POI to better address this condition. Accumulating evidence suggests that mitochondrial events, which contribute to negative effects such as changes in the mitochondrial metabolic pathway, the triggering of inflammatory responses, the disruption of reactive oxygen species (ROS) homeostasis, and the induction of cell apoptosis, can eventually lead to the further development of POI [[41]6–[42]8]. Additionally, dysfunctional mitochondria are implicated in driving ovarian aging through mechanisms such as impaired mitochondrial dynamics, the accumulation of mitochondrial DNA (mtDNA) mutations, and defective electron transport chain function [[43]9]. Specifically, mitochondria play a crucial role in meeting the energy needs of oogenesis and follicle maturation but are also essential in follicular atresia [[44]10, [45]11]. The deletion of the mitochondrial fusion gene Mfn2 in oocytes leads to impaired oocyte maturation and follicle development [[46]12]. Similarly, targeted deletion of the fission gene Drp1 results in reduced oocyte quality [[47]13]. Further research and exploration of mitochondrial-targeting therapies will offer new insights into the management of POI. The progression of POI is often linked to immune dysregulation [[48]14]. Immune factors affect approximately 10–40% of patients with POI [[49]15]. Abnormalities in the inflammatory response have been associated with impairments in ovarian follicular dynamics. For instance, the enhanced infiltration of T helper 1 (Th1) cells leading to follicle atresia and ovarian dysfunction can be effectively counteracted by regulatory T (Treg) cells [[50]16]. A heightened Th1/Treg ratio and changes in the inflammatory cytokine profile are crucial pathogenic factors in POI [[51]16]. Similarly, an increased CD4 + /CD8 + ratio and decreased abundance of CD8 + /CD57 + T cells have been observed in POI patients, possibly as a result of chronic hypoestrogenism [[52]17]. Furthermore, the presence of ovarian autoantibodies like anti-actinin, anti-nuclear factor, and anti-zona pellucida antibodies, along with alterations in lymphocyte subsets, are directly linked to POI [[53]18, [54]19]. Nonetheless, the specific mechanism underlying immune-related POI remains relatively poorly understood. Research on mitochondrial disorders and immunity has garnered significant interest. Mitochondria play crucial roles in the proinflammatory response, influencing the activation and differentiation of immune cells [[55]20, [56]21]. Moreover, several signals derived from mitochondria are involved in the activation of NLRP3, NLRP6, and other inflammasomes [[57]22]. Studies on immunometabolism have indicated that mitochondrial dysfunction regulates T cell function, contributing to T cell-mediated autoimmune diseases [[58]23]. Additionally, disturbances in mitochondrial metabolism and dynamics affect the immune response by regulating the metabolic activity and function of immune cells [[59]24]. However, the intricate interplay between mitochondrial dysfunction and the immune microenvironment in the context of POI remains poorly understood, with limited studies addressing this dynamic relationship. In our study, we investigated how mitochondrial-related genes contribute to the development of POI and their correlation with the immune microenvironment. We utilized data from the NCBI Gene Expression Omnibus (GEO) database, specifically [60]GSE128240, [61]GSE233743, and [62]GSE39501, along with the MitoCarta 3.0 database. Additionally, we examined gene expression signatures using single-cell RNA sequencing data from the [63]GSE200612 dataset. Furthermore, we conducted an analysis to identify potential drug candidates for POI based on the cMap database. Our objectives were to identify novel biomarkers, provide fresh perspectives on therapeutic strategies, and improve the management of POI. By integrating these datasets, we aimed to elucidate the pathogenesis of POI, understand the immune microenvironment, and potentially identify effective treatment options for this condition. Methods Data collection RNA sequencing (RNA-seq) data (accession numbers [64]GSE128240 and [65]GSE233743) and microarray data (accession number [66]GSE39501) were obtained from the GEO database ([67]http://www.ncbi.nlm.nih.gov/geo/). A total of 29 samples, including 14 POI samples and 15 normal tissue samples, were included in the analysis. The [68]GSE128240 dataset contains ovarian tissue from 3 cyclophosphamide (CTX)-treated mice and 3 control mice. The samples of the [69]GSE233743 dataset were obtained from the ovarian tissues of 8 POI model mice and 9 control mice. The [70]GSE39501 dataset contains ovarian tissue from 3 POI model mice and 3 WT control mice. The expression values were preprocessed using the R package “limma”. The fastq raw data was downloaded from the GEO database, first cleaned it with fastp to get clean data, then re-compared it to the standard genome of mice with hisat2, then annotated and quantified the genes with string tie to get the counts and FPKM values of each gene, and finally screened the counts values for differences with deseq2. The single-cell RNA sequencing (scRNA-seq) data of the mice samples (accession number [71]GSE200612) were generated using 10 × Genomics. For this study, we analyzed data from 8-week-old and 12-month-old mice. The specific details of each dataset are presented in Additional file 1: Table S1. Identification of differentially expressed genes (DEGs) The R package “limma” (version 3.5.1) was used to identify DEGs between normal and POI ovarian tissues, and all identified DEGs met p < 0.05 and |logFC|> 1. The resulting DEGs were visualized by Volcano Plot using the R package “ggplot2 ([72]https://www.ncbi.nlm.nih.gov/geo/geo2r/)” and Heatmaps using the R package “pheatmap”. Functional enrichment analysis To reveal the characteristics of the DEGs, Gene Set Enrichment Analysis (GSEA), Gene Ontology (GO) analysis, and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis, were performed using the “clusterProfiler” package. p < 0.05 and p.adj < 0.05 were used as cut-off values to determine the significantly enriched pathways. The results were visualized with ggplot2. Principal component analysis (PCA) Three POI datasets were analyzed using the R package “psych” to obtain the whole clustering profile. Highly variable (variation coefficient > 1) genes were included for PCA and the results were visualized with ggplot2. Identification of mitochondrial-related DEGs (MitoDEGs) A total of 1140 mitochondrial-related genes were obtained from the MitoCarta 3.0 database ([73]http://www.broadinstitute.org/mitocarta). MitoDEGs were screened by crossing the DEGs of interest from each dataset with mitochondrial-related genes, and the results were visualized via Venn diagrams using the R package “VennDiagram” and Heatmaps using the R package “pheatmap.” Analysis of protein‒protein interactions (PPI) PPI networks were constructed using the STRING database ([74]https://string-db.org/) based on the MitoDEGs. The hub-MitoDEGs were selected via the plug-ins CytoHubba ([75]https://apps.cytoscape.org/apps/cytohubba) and MCODE ([76]https://apps.cytoscape.org/apps/mcode) as implemented by Cytoscape software (version 3.8.1). Specifically, for the protein interaction subnetworks screened via the MCODE plugin, the parameter settings of the screening process included Degree Cutoff: 2, Max Depth: 100, K-Core: 2, and Node Score Cutoff: 0.2. Subsequently, the CytoHubba plugin was applied to select hub genes within the PPI network with a Matthews correlation coefficient ≥ 60. By combining the results, the top ten hub-MitoDEGs were selected. Bioinformatic evaluation of mitochondrial characteristics in POI Pearson correlation analysis was applied to assess the interrelationship between the hub-MitoDEGs and mitochondrial respiratory chain complex genes [[77]25], mitophagy genes (GOBP_MITOPHAGY, gsea-msigdb.org), mitochondrial dynamics-related genes (MitoCarta3.0), and mitochondrial metabolism-related genes (MitoCarta3.0). The results were visualized by Heatmaps using the R package “pheatmap”. Exploration of immune characteristics in POI The Rank-In ([78]http://www.badd-cao.net/rank-in/submission.html) was applied in the standardized preprocessing of the gene matrices of the [79]GSE128240, [80]GSE233743, and [81]GSE39501 original datasets. Further immune infiltration analysis was explored based on the normalized gene expression matrix. The formatted data were imported to the Immune Cell Abundance Identifier for mouse (ImmuCellAI-mouse, [82]https://guolab.wchscu.cn/ImmuCellAI-mouse/), which was designed to estimate the infiltration level of 36 immune cells on the basis of the mouse transcriptome data generated by RNA-seq or microarray. The main process includes the following three steps, first calculate the expression deviation matrix, then calculate the single sample Gene Set Enrichment Analysis enrichment score in each cell layer separately. Subsequently, normalize the calculated enrichment scores to the final abundance of immune cells [[83]26]. The Wilcoxon rank sum test was used for group comparisons. The correlations between the hub-MitoDEGs/mitochondrial-related genes and immunocytes, immune checkpoints [[84]27], and immune-related genes (ImmPort, [85]https://www.immport.org/home) were analyzed through Pearson correlation analysis. Clinical samples Granulosa cells (GCs) from POI patients (POI) and healthy individuals (healthy controls, HC) were obtained from women undergoing in vitro fertilization or intracytoplasmic sperm injection and embryo transfer at the Reproductive Medical Center of the Fourth Affiliated Hospital of Jiangsu University (Zhenjiang Maternity and Child Health Care Hospital) from June 2022 to October 2023. The inclusion criteria for the POI group in this study were patients who met the European Society of Human Reproduction and Embryology guideline, including (1) age < 40 years, (2) basal follicle-stimulating hormone (FSH) level > 25 IU/L on two occasions > 4 weeks apart, and (3) oligo/amenorrhea > 4 months [[86]28]. The clinical characteristics of all participants are listed in Additional file 1: Table S2. Establishment of an animal model of POI All animal experiments were performed according to the Institutional Animal Care and Use Committee (IACUC) of the Fourth Affiliated Hospital of Jiangsu University and were approved by the laboratory guidelines for the ethical review of animal welfare. Six-week-old C57BL/6 J female mice were obtained from Jiangsu University Experimental Animal Center. The C57BL/6 J mice were intraperitoneally injected with 50 mg/kg CTX for 14 consecutive days to establish the POI model. Western blot RIPA peptide lysis buffer (Solarbio, China) supplemented with 1% protease inhibitors (Cell Signaling Technology, USA) was applied to extract the protein of GCs and POI mouse ovaries. Proteins were separated by SDS‒PAGE and then transferred to PVDF membranes (Merck KGaA, Germany). Then the PVDF membranes were blocked using 5% skim milk for 2 h. Subsequently, after immunoblotting with primary antibodies against CPT1A (1:1000, BS7744, Bioworld, China) and GAPDH (1:10,000, 10,494–1-AP, Proteintech, China) overnight at 4 °C, the sections were washed in TBST and then incubated with a horseradish peroxidase-conjugated secondary antibody (1:10,000, BL004A, Biosharp, China) for 2 h at room temperature. Finally, the results were visualized via a SuperPico ECL Master Mix (Vazyme, China). RT‒qPCR analysis Ovarian tissues were collected from each group. Total RNA was extracted using TRIzol reagent (Thermo Fisher, USA) according to the manufacturer’s instructions. For mRNA quantification, RNA was converted to cDNA using a reverse transcription kit (Vazyme, China). RT‒qPCR was performed using ChamQ Universal SYBR qPCR Master Mix (Vazyme, China). The primers used are listed in Additional file 1: Table S3. Identification of potential drug candidates Potential therapeutic drugs for POI were identified using the Connectivity Map (cMap, [87]https://clue.io/), which uses an algorithm called the weighted connectivity map score that allows for the calculation of potential gene expression profiles for each agent when applied to specific cell lines. The screened hub genes were input into the cMap database as index gene profiles, the 3D structures were downloaded from the National Institutes of Health database ([88]https://pubchem.ncbi.nlm.nih.gov/compound), then the analysis running automatically. Positive-scoring drug imply high similarity to the input genes, while negative-scoring drugs have opposite or dissimilar roles from input genes and are expected to be candidates for disease management [[89]29, [90]30], the top 9 highly negative-scoring drugs were selected for our study. To explore the specific mechanism of these drug candidates, molecular docking simulations were conducted with AutoDock Vina software (version 1.1.2), and the results were visualized with PyMOL software (version 2.5.7). Statistical analysis Quantitative data are presented as the mean ± standard deviation. Statistical analysis was performed using GraphPad Prism version 9 (GraphPad Software, San Diego, CA, USA). Student’s t test was used to analyze the significant differences (p values) between any two groups. All experiments were conducted at least in triplicate, and the significance level was set at p < 0.05. Results Differentially expressed genes associated with POI and functional enrichment analysis The flowcharts of this study were presented in Fig. [91]1 and Additional file 1: Fig. S1. Three publicly available datasets, namely, [92]GSE39501, [93]GSE128240, and [94]GSE233743, were used for analysis. PCA revealed good separation between the control and POI group (Additional file 1: Fig. S2a). Subsequently, the genes exhibiting significant changes in expression were visualized as heatmaps and volcano plots. A total of 596 differentially expressed genes (DEGs), including 360 upregulated genes and 236 downregulated genes, were detected in [95]GSE39501; 2378 DEGs, including 1343 upregulated genes and 1035 downregulated genes, were detected in [96]GSE128240; and 4714 DEGs, including 1171 upregulated genes and 3543 downregulated genes, were detected in [97]GSE233743 (3484 DEGs in Exp. 1, 755 upregulated and 2729 downregulated; 1971 DEGs in Exp. 2 and 3, 501 upregulated and 1470 downregulated) (Fig. [98]2a–h). Fig. 1. [99]Fig. 1 [100]Open in a new tab Flow diagram of the study Fig. 2. [101]Fig. 2 [102]Open in a new tab DEGs and enrichment analysis of POI. a–d Volcano plot of DEGs between POI and control (CON) groups in [103]GSE39501 (a), GSE 128240 (b), and [104]GSE233743(c, d), red: up-regulated, blue: down-regulated. e–h Heatmap of DEGs between POI and control (CON) groups in [105]GSE39501 (e), GSE 128240 (f), and [106]GSE233743(g, h). i–j GSEA analysis for [107]GSE39501, the top 5 GO terms (i) and KEGG pathways (j) based on the normalized enrichment score (NES) between WT and POI groups. k GSEA analysis for GSE1282401, the top 5 GO terms based on the NES between WT and POI groups GSEA was applied to explore significant functional terms between POI and control samples in each dataset. Results suggested that the DEGs mainly participated in GO terms related to cell cycle, cell division, phosphorylation, positive regulation of transcription by RNA polymerase II, positive regulation of transcription, DNA − templated, and enriched in pathways including FoxO signaling pathway, oocyte meiosis, progesterone − mediated oocyte maturation, regulation of actin cytoskeleton, signaling pathways regulating pluripotency of stem cells in the [108]GSE39501 dataset; the DEGs mainly participated in GO terms related to anterior/posterior pattern specification, glucose homeostasis, proximal/distal pattern formation, skeletal system development, sterol biosynthetic process, and enriched in pathways including adipocytokine signaling pathway, bile secretion, hepatocellular carcinoma, steroid biosynthesis, Wnt signaling pathway in the [109]GSE128240 dataset; the DEGs mainly participated in GO terms related to cholesterol biosynthetic process, cholesterol metabolic process, steroid biosynthetic process, steroid metabolic process, sterol biosynthetic process, and enriched in pathways including axon guidance, DNA replication, fatty acid metabolism, steroid biosynthesis, terpenoid backbone biosynthesis in the [110]GSE233743 Exp. 1 dataset; the DEGs mainly participated in GO terms related to apoptotic process, axon guidance, cellular response to calcium ion, fear response, positive regulation of endothelial cell proliferation, and enriched in pathways including axon guidance, fatty acid metabolism, IL − 17 signaling pathway, NF − kappa B signaling pathway, osteoclast differentiation in the [111]GSE233743 Exp. 2 and 3 dataset (Fig. [112]2i–k and Additional file 1: Fig. S2b–f). Moreover, the DEGs predominantly participated in immune system processes, the cell cycle, multicellular organism development, lipid metabolic processes, and ion transport, as indicated by the GO analysis; additionally, they were enriched in metabolic pathways, signaling pathways regulating pluripotency of stem cells, and neuroactive ligand − receptor interactions according to KEGG enrichment analysis (Additional file 1: Fig. S3–4). These findings provide insight into the molecular mechanisms underlying POI and highlight key pathways and processes associated with this condition, potentially paving the way for the identification of novel biomarkers and therapeutic targets. Mitochondrial-related DEGs in POI Next, we determined the mitochondrial-related DEGs (MitoDEGs) based on the MitoCarta 3.0 database. A total of 119 MitoDEGs (33 upregulated and 86 downregulated, genes with inconsistent expression in various datasets were removed) were filtered. Specifically, 25 MitoDEGs were screened from the [113]GSE39501 cohort, 9 upregulated and 16 downregulated; 73 MitoDEGs were screened from the [114]GSE128240 cohort, 18 upregulated and 55 downregulated; and 93 MitoDEGs were screened from the [115]GSE233743 cohort, 40 upregulated and 53 downregulated (82 MitoDEGs in the [116]GSE233743 Exp. 1, 35 upregulated and 47 downregulated, and 23 MitoDEGs in the [117]GSE233743 Exp. 2 and 3, 8 upregulated and 15 downregulated) (Fig. [118]3a–d). The overlap between MitoDEGs is represented via Venn diagrams (Additional file 1: Fig. S5). Fig. 3. [119]Fig. 3 [120]Open in a new tab MitoDEGs and enrichment analysis of POI. a–d Clustered heatmap of DEGs both in MitoCarta3.0 and [121]GSE39501 (a), GSE 128240 (b), and [122]GSE233743 (c, d). e–f GO analysis of down-regulated MitoDEGs (e) and up-regulated MitoDEGs (f). g–h KEGG pathway enrichment analysis of down-regulated MitoDEGs (g) and up-regulated MitoDEGs (h) GO enrichment analysis indicated that the most relevant GO terms associated with the upregulated MitoDEGs included triglyceride metabolic process and mitochondrial respiratory chain supercomplex assembly, while the most relevant fatty acid metabolic process and lipid metabolic process terms were enriched for the downregulated MitoDEGs (Fig. [123]3e–f). The metabolic pathways were consistently notable, as confirmed by KEGG analysis (Fig. [124]3g–h). These findings suggest significant alterations in mitochondrial-related processes in POI, providing insights into the potential mechanisms underlying this condition and suggesting new avenues for therapeutic interventions and biomarker development. PPI network construction and hub-MitoDEG identification A PPI network of the 119 MitoDEGs was constructed using the STRING database (Fig. [125]4a). The MCODE plugin was applied to identify gene cluster modules, and the CytoHubba plugin was used to score the selected genes (Fig. [126]4b–c). Ten candidate hub-MitoDEGs, namely, Hadhb, Cpt1a, Cpt2, Mrpl12, Mrpl51, Mrps7, Eci1, Acaa1b, Acacb, and Acss1, were identified. To further verify our findings, a POI mouse model was constructed, and human GCs were collected. RT‒qPCR suggested that the expression of Hadhb and Cpt1a increased while that of Mrpl12, Mrps7, and Acacb decreased in both the POI model and GCs (Fig. [127]4d and f). Given that Cpt1a showed the most significant differences, western blot analysis was performed, and the results suggested that the Cpt1a protein level was increased in both POI mouse ovaries and human GCs compared with that in the control group (Fig. [128]4e and g). Taken together, our results support the relevance of these genes in the pathophysiology of POI. Fig. 4. [129]Fig. 4 [130]Open in a new tab Construction and validation of hub-MitoDEGs. a PPI network of MitoDEGs. b A key module with 15 genes as the key gene was screened by MCODE. c Top 25 key genes featured by CytoHubba. d RT-qPCR validation of ten hub-MitoDEGs in mice ovaries (n = 3). e Western blot analysis of Cpt1a protein in mice ovaries. f RT-qPCR validation of ten hub-MitoDEGs in GCs from POI patients (n = 50) and control group (n = 50). g Western blot analysis of CPT1A protein in GCs from POI patients and control group. (****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05) Exploration of the specific expression status of hub-MitoDEGs scRNA-seq data from 8-week-old and 12-month-old C57BL/6 J mice were obtained. The major cell types of the ovary were annotated via the nonlinear dimensionality reduction algorithm tSEN. Briefly, we identified ten clusters of cells and distinguished seven major cell types: GCs (three clusters: 0, 2, and 8), stromal cells (Cluster 1), endothelial cells (Cluster 3), immune cells (two clusters: 4 and 9), smooth muscle cells (Cluster 5), epithelial cells (Cluster 6), and theca cells (Cluster 7) (Fig. [131]5a). Except for GCs, a greater cell proportion was found in the aging group (Fig. [132]5b). Among the upregulated genes, Cpt1a was predominantly expressed in the GCs of the aging group, while Hadhb was mostly expressed in theca cells, corresponding to Clusters 2 and 7. In contrast, the downregulated genes Mrpl12, Mrpl51, and Mrps7 were expressed mainly in GCs, followed by in theca cells and in the epithelium of the control group, while Eci was expressed mostly in theca cells. However, the expression pattern of Cpt2 was not consistent with our previous results, and the expression of Acaa1b, Acacb, and Acss1 could not be detected in this dataset (Fig. [133]5c–i, Additional file 1: Fig. S6). Our results provide insights into the specific expression profiles of hub-MitoDEGs across distinct cell populations within the ovary, shedding light on potential cellular functions and interactions in the context of aging and ovarian health. Fig. 5. [134]Fig. 5 [135]Open in a new tab Cell-specific expression of hub-MitoDEGs. a tSNE cluster map revealing 10 specific clusters representing the major ovarian cell types. b The proportion of different ovarian cell types. c Violin plots showing expression of top 10 hub-MitoDEGs for each cluster. d–i tSNE cluster map showing expression of hub-MitoDEGs. Red lines label the main clusters of interest Mitochondrial characteristics in POI Spearman correlation coefficients between respiratory chain complex genes and hub-MitoDEGs were computed. The results revealed positive correlations of respiratory chain-related genes with Cpt1a, Mrpl12, Mrpl51, Mrps7, and Eci1 expression and negative correlations with Hadhb and Cpt2 expression (Fig. [136]6a–e). Fig. 6. [137]Fig. 6 [138]Open in a new tab Characterization of the interaction between mitochondrial function and hub-MitoDEGs. a–e Correlation between hub-MitoDEGs and mitochondrial respiratory chain complex genes. f–g Correlation between hub-MitoDEGs and mitochondrial fission- (f) and fusion-related genes (g). h Correlation between hub-MitoDEGs and mitophagy-related genes. (****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05) We subsequently investigated the potential metabolic features of mitochondria (Fig. [139]7, Additional file 1: Fig. S7). Seven hub-MitoDEGs (Hadhb, Cpt1a, Cpt2, Mrpl12, Mrpl51, Mrps7, and Eci1) were predominantly associated with genes involved in lipid, amino acid, carbohydrate, metal and cofactor, nucleotide, and vitamin metabolism. These findings are in line with the functional enrichment results mentioned above. Additionally, mitochondrial fission/fusion- and mitophagy-related genes mainly interacted with the seven hub-MitoDEGs mentioned above (Fig. [140]6f–h). These discoveries offer a significant understanding of the metabolic attributes and functional links of crucial mitochondrial genes in POI. Fig. 7. [141]Fig. 7 [142]Open in a new tab Correlations between hub-MitoDEGs and mitochondrial metabolism. a Correlation between hub-MitoDEGs and lipid metabolism-related genes. b Correlation between hub-MitoDEGs and amino acid metabolism-related genes. c Correlation between hub-MitoDEGs and carbohydrate metabolism-related genes. (****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05) Immune cell infiltration analysis in POI Thirty-six immune cell infiltrates were analyzed in the control and POI groups by the ImmuCellAI algorithm. Among them, 17 cell types exhibited significant differences. Monocytes and germinal center B were both less abundant in the POI group, while multiple macrophages (including macrophages, M1 macrophages, and M2 macrophages) and T cells (including T cells, CD4 T cells, CD8 T cells, CD4 Tm cells, T helper cells, CD8 Tc cells, and CD8 Tem cells) were more abundant (Fig. [143]8a–b). Multiple correlations between the infiltrating immunocytes were subsequently observed. According to the immune scores, the strongest synergistic effects were between CD8 T cells and CD8 Tc cells, CD8 T cells and CD8 Tcm cells, macrophages and M1 macrophages, and macrophages and M2 macrophages. However, the competitive effects between CD8 Tcm cells and granulocytes and between M1 macrophages and monocytes were the strongest, followed by those between macrophages and monocytes and between CD8 Tc cells and granulocytes (Fig. [144]8c). Fig. 8. [145]Fig. 8 [146]Open in a new tab Characterization of immune infiltration in POI. a The boxplot of the immune cell proportions. b Stacked bar chart of the immune cell. c The correlation matrix of immune cell proportions. d Correlation between hub-MitoDEGs and 36 immune cells. (****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05) Correlations between the hub-MitoDEGs and immune status Subsequently, the potential relationships between the hub-MitoDEGs/mitochondrial-related genes and immunocytes were assessed. Specifically, Hadhb was positively correlated with the infiltration of T helper cells, CD8 T cells, CD8 Tcm cells, and macrophages but was negatively associated with the infiltration of granulocytes. Cpt1a was positively correlated with the infiltration of basophils and follicular B cells but was negatively associated with NKT, pDC, MoDC, and dendritic cells. Cpt2 was positively correlated with the infiltration of naïve CD4 T cells, CD4 T cells, and T helper cells but was negatively associated with B cells. Mrpl12 was positively correlated with monocyte infiltration but was negatively associated with pDC infiltration. Mrpl51 was positively correlated with the infiltration of neutrophils and Tregs but negatively associated with the infiltration of dendritic cells. Mrps7 was positively correlated with monocyte infiltration but was negatively associated with M2 macrophages, macrophages, and M1 macrophages. Eci1 was positively correlated with neutrophil infiltration but was negatively associated with dendritic cells. Additionally, Acaa1b, Acss1, and Acacb were positively correlated with the infiltration of marginal zone B cells but were negatively associated with the infiltration of eosinophils and memory B cells (Fig. [147]8d). In addition, mitochondrial respiratory genes were positively correlated with monocyte counts but negatively correlated with the numbers of various macrophages and T cells (Additional file 1: Fig. S8a–e). Strong associations were also observed between mitophagy-related genes and immunocytes (Additional file 1: Fig. S8f). The Spearman method was further used to evaluate the relationships between the hub MitoDEGs and immune checkpoints (Fig. [148]9a) or immune-related genes (Fig. [149]9b–c, Additional file 1: Fig. S9–11). In general, 7 hub genes, except for Acaa1b, Acacb, and Acss1, were significantly correlated with genes involved in immune checkpoints, antimicrobial agents, the BCR signaling pathway, chemokine receptors, chemokines, interleukin receptors, interleukins, cytokines, cytokine receptors, the TCR signaling pathway, and natural killer cell cytotoxicity. Fig. 9. [150]Fig. 9 [151]Open in a new tab Identification of the correlation between immune-related genes and hub-MitoDEGs. a Correlation between hub-MitoDEGs and immune checkpoints. b Correlation between hub-MitoDEGs and antimicrobials-related genes. c Correlation between hub-MitoDEGs and BCR signaling pathway-related genes. (****p < 0.0001, ***p < 0.001, **p < 0.01, *p < 0.05) cMap analysis and molecular docking cMap analysis revealed nine candidate compounds—calyculin, amodiaquine, eudesmic acid, cefotaxime, BX-912, prostratin, SCH-79797, HU-211, and pizotifen (Fig. [152]10a). Molecular docking simulations were employed to elucidate the therapeutic mechanisms of these drugs. Calyculin ranked first in the cMap database, while its 3D structure was unattainable from PubChem. Amodiaquine, a histamine receptor agonist with the second highest score, was selected to dock with the proteins encoded by the hub-MitoDEGs due to its involvement in follicular development and ovulation [[153]31]. The binding energies between amodiaquine and Hadhb, Cpt1a, Cpt2, Mrpl51, Eci1, Acacb, and Acss1 were − 7, − 8.5, − 8, − 6.1, − 7.6, − 6.2, and − 8.2 kcal/mol, respectively; representative docking results can be seen in Fig. [154]10b–h. Fig. 10. [155]Fig. 10 [156]Open in a new tab Prediction of candidate drugs. a Top 9 therapeutic drug candidates generated by cMap. b–h The molecular docking simulations of amodiaquine to Hadhb (b), Cpt1a (c), Cpt2 (d), Mrpl51 (e), Eci1(f), Acacb (g), and Acss1 (h) Discussion Infertility affects approximately 186 million individuals globally [[157]32]. POI is a significant contributor to female infertility and has garnered increasing attention [[158]33]. Understanding the intricate mechanisms that underlie POI is crucial for early detection, prevention, and personalized treatment. Because available expression arrays for POI patients are limited, we acquired and analyzed three POI-related GEO datasets from mice in this study. Here, we focused on identifying 119 MitoDEGs in POI and screening ten key hub MitoDEGs. Previous research supports the influential role of the mitochondria in regulating oogenesis, follicle maturation, ovulation, and embryogenesis [[159]34]. Dysfunctional mitochondria and disrupted energy metabolism can lead to granulosa cell apoptosis, follicular atresia, and, subsequently, POI development [[160]8]. Furthermore, maintaining a normal immune environment is essential for ovarian development, and disturbances in immune cells, particularly macrophages and dendritic cells, along with imbalances in inflammatory cytokines, are closely associated with POI progression [[161]15]. Variations such as reduced effector suppressor/cytotoxic lymphocyte counts, elevated NK cell numbers and activity, and augmented production of autoantibodies by B cells have been linked to autoimmune ovarian damage [[162]35]. Our aim was to elucidate the pathobiology of mitochondrial function, the immune characteristics in individuals with ovarian impairment, and their interplay in POI development and, consequently, to enhance the clinical effectiveness of POI treatment. In the initial step, we collected human GCs and developed a POI model to validate our findings. The expression patterns of Hadhb, Cpt1a, Mrpl12, Mrps7, and Acacb were consistent with our previous bioinformatic analysis. Given that GCs are known to regulate follicle initiation and development, any dysfunction in these cells may represent a crucial pathological aspect of POI [[163]36]. To further explore this phenomenon, we examined the cell-specific expression profiles of the hub-MitoDEGs in conjunction with the scRNA-seq data. Our results revealed a focused distribution of the hub-MitoDEGs, with upregulated genes such as Hadhb and Cpt1a and downregulated genes such as Mrpl12, Mrpl51, Mrps7, and Eci1 predominantly present in GCs and theca cells. Interestingly, another bioinformatics study conducted using peripheral blood mononuclear cells from POI patients highlighted that the inflammatory response primarily impacts steroid-producing cells, particularly the theca layers and GCs within preovulatory follicles [[164]37]. These results aligned with our own results, suggesting that GCs and theca cells could be the primary targets affected in ovarian dysfunction. Our research validated seven hub genes (Hadhb, Cpt1a, Cpt2, Mrpl12, Mrpl51, Mrps7, and Eci1) that exhibited significant correlations with mitochondrial function or immunity, underscoring the crucial involvement of irregular mitochondrial metabolism in the pathology of POI. Mitochondrial β-oxidation represents a major pathway for ATP production and is vital for oocyte maturation [[165]38]. Key enzymes involved in β-oxidation, such as Hadhb, Cpt1a, Cpt2, and Eci1, oversee the regulation of fatty acid metabolism [[166]39–[167]41]. Any abnormalities in these genes could compromise ovarian function and lipid metabolism. Notably, a decrease in Cpt1a within placenta-derived mesenchymal stem cells has been linked to alterations in energy metabolism and the reversal of senescence [[168]42]. Studies on liver injury have suggested that Cpt1a might exacerbate ROS-induced oxidative stress and inflammation [[169]43]; hence, silencing Cpt1a may protect against oxidative injury. Furthermore, mitochondrial ribosomes, which are responsible for translating subunits of the oxidative phosphorylation (OXPHOS) machinery, play critical roles in ovarian development and function [[170]44]. Mrpl12 contributes to maintaining mitochondrial homeostasis by regulating mitochondrial biosynthesis and activating OXPHOS [[171]45], while Mrpl51 is crucial for maintaining mitochondrial redox balance, mtDNA stability, and integrity [[172]46]. Additionally, Mrps7 has been identified as a diagnostic gene for syndromic POI [[173]47]. Furthermore, previous studies have provided preliminary insights into the molecular regulatory mechanisms of these hub genes. Specifically, estrogen receptor α, a transcription factor, has been reported to interact with Hadhb, therefore contributing to the regulation of lipid metabolism [[174]40]. CircHIPK3/miR-124-3p targets Cpt1a-driven fatty acid oxidation (FAO) to regulate angiogenesis in early-onset preeclampsia [[175]48]. MiR-26, which inhibits fatty acid oxidation and synthesis and decreases the expression of Cpt1a and Acacb, has currently been identified as a promising target for atherosclerosis therapy [[176]49]. In addition, the transcription factors E2F2 and E2F1 both target and inhibit Cpt2 expression, thereby repressing FAO and promoting the formation of a high-lipid environment essential for hepatocarcinogenesis [[177]50]. A recent study suggested that the PI3K/mTOR/Yin Yang 1 pathway coregulates Mrpl12 expression, which is involved in the malignant progression of hepatocellular carcinoma [[178]45]. However, the biological functions of these hub genes in POI have yet to be explored, and future studies will focus on the construction and validation of interested-hubgenes-transcription factors-microRNAs regulatory network. Together, these findings support the potential of these genes as biomarkers for early detection and timely intervention in POI. However, further investigations are required to elucidate the specific functions of these hub genes in the context of POI. Subsequently, we explored the correlation between hub biomarkers and mitochondrial-related genes. ATP production is facilitated by mitochondrial electron transport complexes [[179]51]. Our observations revealed a positive correlation of respiratory chain-related genes with Cpt1a, Mrpl12, Mrpl51, Mrps7, and Eci1 but a negative correlation with Hadhb and Cpt2. Furthermore, altered mitochondrial dynamics are recognized as a potential mechanism contributing to compromised mitochondrial function [[180]52]. Significantly, seven hub-MitoDEGs (Hadhb, Cpt1a, Cpt2, Mrpl12, Mrpl51, Mrps7, and Eci1) were closely associated with mitochondrial dynamics, highlighting their role in coordinating mitochondrial fission and fusion. Mitophagy, which is responsible for selectively degrading damaged mitochondria to uphold immune responses, was the focal point of our study [[181]53]. We observed strong correlations between mitophagy-related genes and the aforementioned seven hub-MitoDEGs and immune cells, particularly interactions with macrophages, monocytes, T cells, CD4 T cells, and CD8 T cells. Besides, metabolomics studies have shed light on how metabolic dysfunction might underpin the long-term complications of POI, notably, glucose and lipid metabolism disorders have been closely linked to POI [[182]8], which aligns well with our findings from the GO, KEGG, and GSEA analysis. Collectively, these findings underscore the importance of safeguarding mitochondrial function during interventions targeting POI, suggesting significant benefits for various metabolic pathways and immune responses. Moreover, our research further evaluated the immune characteristics of POI, revealing distinct immune patterns between the POI group and the control group. Notably, we observed a greater abundance of cDC2s, multiple macrophages, and various T cells in the POI group than in the control group, while monocyte and germinal center B-cell numbers were notably reduced. In line with our results, monocytes were found to be exhausted in POI, suggesting a heightened innate immune response and compromised classical monocyte functionality [[183]37]. Monocytes, characterized as the most dynamic cell population in the mononuclear phagocyte system, are polarized to the classically activated macrophage (M1-like) phenotype, which induces a persistent immune response, or to alternatively activated macrophage (M2-like) repair macrophages, which are related to tissue homeostasis maintenance and inflammation resolution [[184]54–[185]56]. Moreover, in solid tumors, data from several studies revealed that monocytes are recruited and then differentiate into tumor-related macrophages, which accelerate the malignant progression and angiogenesis of tumors [[186]57]. Of particular importance is the role of macrophages as crucial regulatory and effector cells in immune responses, and their activation has been demonstrated to induce ovarian failure in zebrafish [[187]58, [188]59]. Furthermore, the phenotypic transition between M1-like and M2-like macrophages, facilitated by Hadhb-mediated fatty acid oxidation, presents a promising avenue for therapeutic interventions in inflammatory diseases [[189]60]. Mitochondrial function and immunity are inextricably interrelated. A strong linkage between immunocytes and the hub-MitoDEGs was confirmed. In fact, immunomodulatory therapies have previously been reported to normalize the ovarian function of POI patients [[190]61], and the exploration of immune checkpoints is beneficial for improving immunotherapy efficacy [[191]62]. Here, our results revealed interactions among diverse immune checkpoint proteins and Hadhb, Cpt1a, Cpt2, Mrpl12, Mrpl51, Mrps7, and Eci1. Furthermore, various immunocytes were positively correlated with Hadhb and Cpt2 but negatively related to Cpt1a, Mrpl12, Mrpl51, Mrps7, and Eci1. Besides, mitochondria are critical for the regulation of energy production and redox signaling in monocytes [[192]63]. Consistently, our results showed a positive correlation between monocyte enrichment and mitochondrial complex gene expression. These findings suggest that mitochondrial damage may be a potential mechanism impacting the fate and function of immunocytes, which are vital for POI pathogenesis. Finally, using drug information from the cMap database, we identified nine candidate drugs that may improve POI. In detail, the protein phosphatase inhibitor calyculin has been reported to prevent apoptosis by preventing Bax translocation and cytochrome c release [[193]64]. Moreover, transient exposure of oocytes to calyculin-A permits subsequent development and fertilization [[194]65]. The histamine receptor agonist amodiaquine has been routinely prescribed as an antimalarial drug [[195]66]. Interestingly, recent studies have suggested that the application of amodiaquine in diabetic tubulopathy ameliorates mitochondrial disorders by attenuating mitochondrial Nox4 [[196]67]. In addition, histamine and its agonist can induce ovulation in isolated perfused rat ovaries [[197]68]. The bacterial cell wall synthesis inhibitor cefotaxime, a third-generation cephalosporin, has been applied in multiple infections [[198]69]. The pyruvate dehydrogenase kinase (PKD) inhibitor BX-912 induces cell apoptosis and blocks the G[2]/M phase in mantle cell lymphoma cell lines, therefore exerting antiproliferative effects [[199]70]. Another PKD inhibitor, dichloroacetate, has been reported to increase ROS production in M1- and M2- macrophages and induce a metabolic shift from glycolysis to OXPHOS in M1 macrophages [[200]71]. Recent studies have suggested that PKC activators restrict innate immune suppression and promote antigen cross-presentation [[201]72]. Prostratin, a PKC activator, serves as a therapeutic agent that targets HIV viral reservoirs [[202]73]. SCH-79797, a PAR1 antagonist, protects against ischemia/reperfusion injuries in rat hearts and was recently identified as a broad-spectrum antibacterial agent [[203]74, [204]75]. The glutamate receptor antagonist HU-211 has been confirmed to be an effective neuroprotectant [[205]76]. The 5-HT2 receptor antagonist pizotifen is currently applied for vascular headache prevention [[206]77]. However, the application of these drugs in POI treatment is currently lacking; thus, subsequent studies will examine the efficacy of each drug in GCs or in POI animal models. Overall, in the process of integrating and analyzing data related to POI using bioinformatics tools, we gained a deeper understanding of the intricate interplay between mitochondrial function and immune status. Notably, our study identified seven genes (Hadhb, Cpt1a, Cpt2, Mrpl12, Mrpl51, Mrps7, and Eci1) as potential biomarkers for targeted therapy in POI. However, the limitations of our study include the following: (1) clinical POI data from public databases and published studies are currently lacking; hence, future studies should further expand the sample size and perform clinical sample transcriptome sequencing to verify our findings; (2) despite the identification of promising candidate drugs, the lack of basic or clinical trials of these pharmaceuticals indicates the necessity for future in vivo and in vitro experiments; and (3) the specific mechanism of our selected hub biomarkers in mitochondrial function, immunity and their involvement in POI development needs further investigation in actual POI cases and detailed clinical data. Subsequent studies will focus on these directions. Conclusions In brief, through comprehensive bioinformatics analysis, we identified DEGs between the POI and CON groups and subsequently explored the mitochondrial genes related to POI. Functionally, we identified tight links between seven hub-MitoDEGs and mitochondrial function, as well as between immunocyte infiltration and the immune microenvironment, and demonstrated crosstalk between mitochondrial function and immunity. In addition, nine candidate drugs for POI treatment were predicted to be promising pharmacologic agents for future research. Supplementary Information [207]12916_2024_3675_MOESM1_ESM.pdf^ (2.4MB, pdf) Additional file 1. Figures S1-S11. Fig S1. The data processing and flowchart of this study. Fig S2. Results of PCA and GSEA analysis. Fig S3. GO and KEGG analysis of GSE39501 and GSE128240. Fig S4. GO and KEGG analysis of GSE233743. Fig S5. MitoDEGs in POI. Fig S6. Cell-specific expression of Eci1, Acaa1b, Acss1, Acacb. Fig S7. Correlations between hub-MitoDEGs and mitochondrial metabolism. Fig S8. Correlations between mitochondrial related genes and immunocytes. Fig S9. Correlations between hub-MitoDEGs and immune related genes. Fig S10. Correlations between hub-MitoDEGs and cytokines receptors- and cytokines- related genes. Fig S11. Correlations between hub-MitoDEGs and TCR signalling pathway-, TGFb family member receptors-, TGFb family members -, TNF family member receptors-, and TNF family members- related genes. Table S1-S3. Table S1: Sample information for this study. Table S2. Patient information for this study. Table S3. Sequences of the primers used in this study. Acknowledgements