Identification of hypoxia- and mitochondria-related genes associated with recurrent spontaneous abortion: an integrative bioinformatics analysis and preliminary validation

Data collection and gene set definition

RSA-associated datasets (GSE165004 and GSE26787) were obtained from the GEO database. The GSE165004 dataset (GPL16699) contains 24 endometrial tissue samples from patients with RSA and 24 from normal controls, and was used as the training set. The GSE26787 dataset (GPL570) contains five RSA and five control endometrial tissue samples, and served as the validation set.

Using GSEA and the MitoCarta 3.0 database, we identified 2,030 mitochondria-related genes (MRGs) for further analysis. Additionally, 493 hypoxia-related genes (HRGs) were collected from published literature (Yu et al., 2023).

Identification of differentially expressed genes

To identify DEGs between RSA and control samples in the GSE165004 dataset, differential expression analysis was performed using the limma package (v1.38.0) (Li et al., 2022). Raw expression data were processed using Agilent Feature Extraction Software (v11.0.1.1) to extract hybridization signals. Background correction and quantile normalization were subsequently performed using limma to ensure data quality and comparability. DEGs were identified based on the criteria p < 0.05 and —log2FoldChange—> 0.5 (James et al., 2022). The distribution of DEGs was visualized using a volcano plot (ggplot2, v3.4.1) (Wang et al., 2022a) and a heatmap (ComplexHeatmap, v2.14.0) (Ding, Goldberg & Zhou, 2023). Notably, the GSE165004 and GSE26787 datasets were analyzed independently due to differences in microarray platforms (GPL16699 and GPL570, respectively), so no data integration was performed.

Overlap and functional profiling of candidate genes

Candidate genes potentially related to RSA were identified by intersecting DEGs with HRGs and MRGs. To explore possible biological functions and pathways of these genes, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were conducted using the clusterProfiler package (v4.7.1.3) (Yu et al., 2012), with a significance cutoff of p < 0.05. Additionally, protein-protein interaction (PPI) networks were constructed via the STRING database (interaction score ≥ 0.4) to explore molecular interactions among candidate genes (Deng et al., 2023).

Selection of gene signatures via machine learning

After obtaining candidate genes, Least Absolute Shrinkage and Selection Operator (LASSO) and Support Vector Machine Recursive Feature Elimination (SVM-RFE) algorithms were applied using the glmnet (version 4.1.4) (Shi et al., 2023) and caret package (version 6.0.93) (Zhang et al., 2019) to select potential signature genes. The intersection of signature genes derived from two machine-learning algorithms was defined as candidate signature genes. The expression of these genes was further explored, and those showing significant differences between RSA and control groups with consistent trends in both GSE165004 and GSE26787 datasets were selected as hypoxia- and mitochondria-related genes (p < 0.05).

Development of a gene-based risk estimation model

To preliminarily evaluate the potential predictive value of hypoxia- and mitochondria-related genes for RSA risk, a nomogram was developed using the rms package (v6.0-1) (Liu et al., 2021), based on their expression levels. Calibration curve and HL test (p > 0.05) were used to assess model performance. Furthermore, performance of nomogram was further evaluated by decision curve analysis (DCA) which was displayed by rmda package (version 1.6).

Enrichment analysis based on selected gene expression

To investigate the biological pathways associated with the expression of hypoxia- and mitochondria-related genes, RSA samples from the GSE165004 dataset were divided into high- and low-expression groups based on the median expression of each gene. Differential gene expression analysis between the two groups was performed using the limma package. All genes were ranked according to log2 fold change values between high and low expression groups. Gene Set Enrichment Analysis (GSEA) was then conducted using the clusterProfiler package (v4.7.1.3), with the ranked gene list as input. The background gene sets used for enrichment were obtained from the Molecular Signatures Database (MSigDB), including GO terms (c5.go.v2023.2.Hs.symbols.gmt) and KEGG pathways (c2.cp.kegg.v2023.1.Hs.symbols.gmt). Enriched terms with —NES—>1 and p-value < 0.05 were considered statistically significant. Tables S1–S2 include the ranked gene lists, and Tables S3–S4 contain details of the background gene sets used.

Immune cell composition and gene expression correlation

The proportions of 22 immune cell types in the GSE165004 dataset were estimated using the CIBERSORT algorithm. Immune cell types with a score of 0 in over 25% of samples were excluded. Differences in immune infiltration between RSA and control groups were assessed using the Wilcoxon test (p < 0.05). Spearman correlation analysis was further applied to examine the relationship between hub genes and differential immune cells (—r—> 0.3, p < 0.05).

Prediction of regulatory associations and related compounds

To explore potential regulatory mechanisms involving the RSA-associated genes identified in this study, transcription factors (TFs) were predicted using the RNAInter dataset, and those with a confidence score >0.1 were retained. In addition, microRNAs (miRNAs) predicted to target these genes were identified using miRDB, TargetScan, and miRWalk. Key miRNAs were obtained by taking the intersection of miRNAs derived from three databases. To screen candidate compounds potentially linked to these RSA-related genes, we predicted compounds via the Enrichr platform using DSigDB, based on the expression profiles of selected genes. Specifically, these genes were input into Enrichr, and compounds with significant enrichment scores were identified based on their associations with gene expression profiles in DSigDB. The predicted drugs were ranked according to their adjusted p-values, and a subset of top-ranked candidates with statistically significant associations was selected for further analysis. Finally, a compound–gene interaction network was visualized using Cytoscape software (v3.10.1) (Wang et al., 2023).

RT-qPCR assessment of selected gene expression

To validate the expression of hypoxia- and mitochondria-related genes in endometrial tissues, specimens were obtained from five individuals with RSA and five matched healthy controls. This study has obtained the approval of the Ethics Committee of Jiangbin Hospital Guangxi Zhuang Autonomous Region (LW2025-003). All patients included in this study provided written informed consent and granted permission for the use of their medical records for the current research.

Total RNA was extracted from 50 mg of frozen tissue using the TRIzol method (Ambion, Austin, TX, USA). Samples were thoroughly homogenized in one mL of TRIzol and incubated on ice for 10 min; we then added 300 µL of chloroform (Chengdu Goodyear Gum Co., Ltd.). After standing at room temperature for another 10 min, the samples were centrifuged at 12,000 × g for 15 min at 4 °C. The upper aqueous phase was carefully collected, mixed with an equal volume of isopropanol (Chengdu Cologne Chemical Co., Ltd.) to precipitate RNA, washed twice with 75% ethanol, and finally dissolved in 20–50 µL of RNase-free water (Xavier).

Tissues were snap-frozen in liquid nitrogen within 15 min of collection and stored at −80 °C until use. Genomic DNA was removed by on-column DNase I treatment (Thermo Fisher Scientific, Waltham, MA, USA), following the manufacturer’s protocol. RNA concentration and purity were assessed using a NanoPhotometer N50 (see Table S5). No-reverse transcriptase (no-RT) controls were included to confirm the absence of genomic DNA contamination.

First-strand cDNA was synthesized using the SweScript First Strand cDNA Synthesis Kit (Xavier) in a 20 µL reaction containing four µL of 5× Reaction Buffer, one µL of Oligo(dT) primer (provided in the kit), one µL of SweScript RT I Enzyme Mix, and 0.1 ng to five µg of total RNA (amount adjusted based on concentration). Reverse transcription was performed under the following conditions: 25 °C for 5 min, 50 °C for 15 min, and 85 °C for 5 s. The resulting cDNA was diluted 5–20-fold with nuclease-free water and stored at −80 °C. To assess potential inhibition, cDNA samples were tested in serial dilutions (5 ×, 10 ×, 20 ×), and no inhibition was observed.

RT-qPCR was performed using the CFX Connect™ Real-Time PCR Detection System (Bio-Rad, Hercules, CA, USA). Each 10 µL reaction consisted of 3 µL of diluted cDNA, 5 µL of 2 × Universal Blue SYBR Green qPCR Master Mix (Xavier), and 0.5 µM each of forward and reverse primers. Primers were designed to amplify all known transcript variants of the target genes. Primer specificity was confirmed using NCBI Primer-BLAST. Detailed sequences, accession numbers, and exon–intron boundary information for all primers are listed in Table S6. Amplification curves were used to evaluate PCR efficiency and relative gene expression levels, while melting curve analysis verified primer specificity (see Fig. S1). No-template control (NTC) reactions were included in each run. No amplification was observed in the NTCs, confirming the absence of contamination or primer-dimer formation. The thermal cycling program included an initial denaturation at 95 °C for 1 min, followed by 40 cycles of 95 °C for 20 s, 55 °C for 20 s, and 72 °C for 30 s, with subsequent melting curve analysis from 65 °C to 95 °C at 0.5 °C increments every 5 s.

Each sample was analyzed in technical triplicate, with 5 biological replicates per group (Control and RSA). Gene expression was quantified using the 2^ΔΔCt method, with GAPDH as the internal reference gene (H-GAPDH① for the Control group and H-GAPDH② for the RSA group). GAPDH was selected based on prior literature and its consistent expression in endometrial tissue across conditions. Ct values above 35 were excluded, and technical replicates with a coefficient of variation exceeding 5% were considered outliers and removed. Statistical analysis was conducted using unpaired t-tests in GraphPad Prism 5, with significance set at p < 0.05.

Statistical analysis

Statistical analyses were performed using R software (version 4.2.2) and GraphPad Prism 5. DEGs were identified using the limma package; immune cell differences were compared using the Wilcoxon rank-sum test, correlations were assessed by Spearman analysis, and RT-qPCR data were analyzed using unpaired t-tests, with p < 0.05 considered statistically significant.

Screening and analysis of 11 RSA-associated candidate genes

A total of 831 DEGs, including 419 upregulated and 412 downregulated genes, were identified between the RSA and control groups in the GSE165004 dataset (Figs. 1A–1B). 11 candidate genes were selected by intersecting the 831 DEGs with 2,030 MRGs and 493 HRGs (Fig. 1C, Table S7). Enrichment analysis of these genes revealed enrichment in 584 GO terms and 16 KEGG pathways. GO terms included response to hypoxia, cellular response to oxygen levels, and response to decreased oxygen, among others (Fig. 1D, Table S8). KEGG pathways associated with these genes included hypertrophic cardiomyopathy, AGE-RAGE signaling in diabetic complications, and HIF-1 signaling, among others (Fig. 1E). In the PPI network, UCP3 showed relatively higher connectivity with other RSA-associated candidate genes (Fig. 1F).

Expression analysis of BSG and TRAK1 as potential RSA-associated genes

To explore genes potentially associated with RSA, machine-learning algorithms were applied. LASSO analysis yielded UCP3, BCL2, BSG, CYBB, TRAK1, and EPAS1 as one set of candidate genes (lambda._min = 0.00924) (Fig. S2A). Additionally, SVM-RFE produced another set of genes—UCP3, BCL2, UCP2, BSG, CYBB, TRAK1, and SLC8A1—corresponding to the point where the model showed relatively higher classification performance (Fig. S2B). By intersecting the potential signature genes from both machine-learning algorithms, five candidate signature genes (UCP3, BCL2, BSG, CYBB, and TRAK1) were selected (Fig. 2A). Further expression analysis showed that BSG and TRAK1 were significantly upregulated in RSA groups in both the GSE165004 and GSE26787 datasets (Figs. 2B–2C). These findings suggest that BSG and TRAK1 are potentially associated with RSA in these cohorts.

Predictive nomogram based on BSG and TRAK1 expression in RSA patients

A nomogram model was constructed using the expression levels of BSG and TRAK1 to explore their potential utility in assessing RSA-related patterns (Fig. 3A). The calibration curve showed visual alignment with the reference line, and the p-value from the HL test was greater than 0.05, suggesting that the model may exhibit a degree of calibration in this dataset (Fig. 3B). DCA results also showed that the nomogram may have potential clinical relevance within this sample context (Fig. 3C).

Enrichment analysis of pathways related to BSG and TRAK1 expression in RSA

GSEA on genes ranked by differential expression between BSG/TRAK1 high- and low-expression groups revealed significant pathway enrichment. For BSG, a total of 1,355 GO terms and 60 KEGG pathways were enriched. Representative GO terms BSG included outer kinetochore, intracellular zinc ion homeostasis, and detoxification of copper ions (Fig. 4A). Enriched KEGG pathways involved DNA replication, mismatch repair, and glycosaminoglycan (GAG) biosynthesis-keratan sulfate (Fig. 4B). In the case of TRAK1, 1,427 GO terms and 38 KEGG pathways were significantly enriched. Key GO terms included serine-type endopeptidase complex, dendritic cell (DC) apoptotic process, and oxalate transport (Fig. 4C). The enriched KEGG pathways included complement and coagulation cascades, natural killer (NK) cell-mediated cytotoxicity, and terpenoid backbone biosynthesis (Fig. 4D).

Immune cell composition and gene correlation patterns in RSA

The proportions of 22 immune cell types across all samples from the GSE165004 dataset are presented in Fig. 5A. Differences between the RSA and control groups were observed in M2 macrophages, monocytes, resting memory CD4+ T cells, CD8+ T cells, gamma delta T cells, and regulatory T cells (Tregs). Specifically, M2 macrophages, CD8+ T cells, and gamma delta T cells (γδ T cells) were more prevalent in the control group, whereas the remaining differentially expressed immune cell types exhibited the opposite pattern (Fig. 5B). Furthermore, a Spearman correlation analysis between the selected genes and differentially infiltrated immune cells revealed a positive correlation between BSG and resting memory CD4+ T cells (r > 0.3, p < 0.05) (Fig. 5C). TRAK1 was negatively correlated with M2 macrophages, CD8+ T cells, and γδ T cells (r <−0.3, p <  0.05), but positively correlated with resting memory CD4+ T cells and Tregs (r > 0.3, p < 0.05) (Fig. 5D).

Gene regulatory associations and predicted compounds related to BSG and TRAK1 in RSA

To explore potential regulatory associations of selected genes in RSA-related datasets, 19 TFs were computationally predicted. Among these, EHMT2 and SUPT5H were predicted to be potentially associated with the regulation of both BSG and TRAK1 (Fig. 6A). Based on the selected genes, 258 miRNAs were retrieved from the miRDB database, 1,043 from TargetScan, and 3,996 from miRWalk, resulting in 122 candidate miRNAs by intersecting the outputs of all three databases. Among these, 11 candidate miRNAs (e.g., hsa-miR-6766-5p, hsa-miR-1237-5p, hsa-miR-3692-3p) were predicted for BSG, and 111 key miRNAs (e.g., hsa-miR-6511b-5p, hsa-miR-921, hsa-miR-1251-3p) for TRAK1 (Fig. 6B). Captopril, valproic acid, and compound 7646-79-9 were computationally associated with BSG and TRAK1 in drug prediction analysis, but no functional link is confirmed (Fig. 6C).

RT-qPCR verification of BSG and TRAK1 expression

RT-qPCR analysis demonstrated a significant increase in BSG expression in the RSA group compared to the control group (p = 0.0010; Fig. 7A). TRAK1 expression showed an increasing trend in the RSA group compared with the control group, though the difference was not statistically significant (p = 0.4415; Fig. 7B).