Gene expression profiles of PE were retrieved from the National Center for Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database (7). Three placental tissue-derived Homo sapiens datasets were analyzed: GSE75010 (GPL6244 platform; 80 PE cases and 77 controls) (8–13), GSE10588 (GPL2986 platform; 17 PE cases and 26 controls) (14), and GSE54618 (GPL10558 platform; 12 PE cases and 12 controls) (15) ( Table 1 ). Probe annotation was performed for GSE75010 and GSE10588 using platform-specific annotation files, followed by dataset merging and batch effect correction via the sva R package (v3.52.0) (16). The combined dataset (97 PE cases and 103 controls) underwent normalization using the R package limma (v3.60.4) (17), with principal component analysis (PCA) (18) confirming effective batch effect removal. The GSE54618 dataset served as an independent validation cohort, processed identically with probe annotation and normalization.

RiboRGs were systematically acquired through a dual-source approach. Firstly, the GeneCards database (19) was queried using the keyword “ribosome biogenesis,” retaining protein-coding genes with a relevance score >5, which yielded 59 candidate genes. Secondly, 331 RiboRGs were retrieved from the published literature using the same keyword on the PubMed website (20). After merging these two gene lists and removing duplicates, a final set of 344 nonrepetitive RiboRGs was generated for subsequent analysis ( Supplementary Table S1 ).

Differentially expressed genes (DEGs) between PE and control groups were identified using the R package limma with thresholds of |log2 fold change (log2FC) | > 0.1 and adjusted p-value (P.adj) < 0.05 (Benjamini-Hochberg correction). Upregulated and downregulated DEGs were defined as log2FC > 0.1 and log2FC < -0.1, respectively, with statistical significance (P.adj < 0.05). Volcano plots were generated using R package ggplot2 (v3.5.1) to visualize differential expression patterns. RiboRDEGs were subsequently identified by intersecting the DEG list with the precompiled 344 RiboRGs, with results visualized through a Venn diagram. Heatmaps generated by R package ComplexHeatmap (v2.20.0) (21) employed Z-score-normalized counts, hierarchical clustering (Euclidean distance, complete linkage), and three-dimensional PCA maps generated by R package rgl (v1.3.1) demonstrated clear separation of PE and controls.

Gene set enrichment analysis (GSEA) (22) was performed on the combined dataset using the R package clusterProfiler (v4.12.6) (23) with the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (24), employing gene set size thresholds of 10–500 genes and significance criteria of P.adj < 0.05 (Benjamini-Hochberg method) and false discovery rate (FDR) < 0.05. To elucidate the biological functions of RiboRDEGs, Gene Ontology (GO) (25), which encompasses biological process (BP), cellular component (CC), and molecular function (MF), and KEGG pathway analysis were performed using clusterProfiler, with the same statistical thresholds. Pathway interaction networks were reconstructed using the R package CBNplot (v1.4.0) (26), where Bayesian networks were inferred through the bnpathplot function by modeling biological pathways as nodes weighted by pathway activity scores derived from RiboRDEG expression profiles.

WGCNA was implemented using the R package WGCNA (v1.73) (27) to identify hub genes among RiboRDEGs. A scale-free topology model was constructed by selecting an optimal soft-thresholding power (β) to maximize network connectivity while minimizing spurious correlations. The adjacency matrix was transformed into a topological overlap matrix (TOM) to quantify gene co-expression similarity, followed by dynamic tree cutting to define gene modules. Module-trait relationships were assessed by calculating Pearson correlation coefficients between module eigengenes (MEs) and PE status, with the most significantly associated module (p < 0.05) selected for downstream analysis. Hub genes were identified as the intersection of RiboRDEGs and genes within the PE-correlated module, visualized via Venn diagrams. Pairwise Spearman correlations among hub genes were computed and displayed using the corrplot R package (v0.94).

A comprehensive machine learning framework comprising (28) 113 prediction models was developed using the combined dataset to identify key feature genes associated with PE risk. Twelve distinct algorithms spanning linear models (Stepglm, Lasso, Ridge, Enet), ensemble methods (XGBoost, RandomForest, GBM), Bayesian approaches (NaiveBayes), hybrid dimensionality reduction & regularization (plsRglm) and supervised learning techniques (SVM, glmBoost, LDA) were implemented through their respective R packages (glmnet, xgboost, randomForest, etc.), with algorithm combinations detailed in Supplementary Table S2 . Model performance was evaluated via 10-fold cross-validation, receiver operating characteristic (ROC) curve analysis, and decision curve analysis (DCA), with diagnostic efficacy quantified by mean area under curve (AUC) values across training and validation cohorts (GSE54618). Models demonstrating AUC >0.9 were prioritized as high-diagnostic-value candidates. Final feature gene selection was guided by consensus across top-performing models, validated through calibration curves and confusion matrices to ensure robustness.

Model interpretability was assessed using SHAP (29) to delineate the contribution of key feature genes to PE risk prediction. SHAP values were computed via the R package kernelshap (v0.7.0), with positive/negative values indicating directional effects on risk (increase/decrease). Global feature importance rankings were derived from mean absolute SHAP values using shapviz (v0.9.5), visualized through bar plots (overall importance) and beeswarm plots (feature value-SHAP value distributions). SHAP interaction values quantified pairwise feature interdependencies, visualized via scatterplots, while waterfall plots generated for representative cases provided localized interpretability of model decisions. This comprehensive SHAP-based interpretability framework quantifies the importance of key feature genes in predictive model decisions in PE.

Key genes were assessed for differential expression between PE and control groups using Mann-Whitney U tests, with results visualized through violin plots. Diagnostic performance was evaluated via ROC curve analysis using the R package pROC (v1.18.5) (30), calculating AUC values to quantify predictive capacity. PPI networks were reconstructed via the GeneMANIA database (31), integrating key genes with functionally associated partners to infer biological modules relevant to PE pathogenesis.

Transcriptional and post-transcriptional regulatory networks were systematically reconstructed to elucidate molecular interactions involving key genes. Transcription factor (TF)-gene interactions were identified using the ChIPBase database (32), retaining TF-gene pairs with combined upstream/downstream supporting samples ≥8. RNA binding protein (RBP) and microRNA (miRNA) interactors were predicted via the ENCORI database (33), applying evidence-based thresholds of clipExpNum ≥10 for RBPs and ≥7 for miRNAs. All interaction networks (TF-gene, RBP-gene, miRNA-gene) were integrated and visualized using Cytoscape (v3.10.2).

Placental immune cell infiltration profiles were quantified via single-sample gene set enrichment analysis (ssGSEA) using the R package GSVA (v1.52.3), with 28 immune cell-specific gene sets derived from established markers (34). Enrichment scores representing relative immune cell abundance were compared between PE and control groups, identifying differentially infiltrated cell types. Spearman correlation matrices generated via the R package linkET (v0.0.7.4) revealed intercellular immune interactions and key gene-immune cell associations (|r| >0.3), visualized through network diagrams. Significant correlations (|r| >0.3) were further validated using scatterplots to delineate linear relationships.

Placental tissue samples were collected from 20 singleton pregnancies (10 PE cases, 10 gestational age-matched controls, each group included five term pregnancy and five preterm pregnancies) undergoing cesarean delivery at the Third Affiliated Hospital of Wenzhou Medical University. PE diagnosis followed ISSHP criteria: sustained hypertension (≥140/90 mmHg) with proteinuria (≥0.3 g/24h) emerging after 20 gestational weeks. Inclusion criteria required maternal age 20–40 years, uncomplicated antenatal course prior to PE onset, and absence of fetal anomalies. Exclusion criteria encompassed pre-existing comorbidities (chronic hypertension, diabetes), acute infections (including COVID-19), gestational diabetes, fetal congenital disorders, and exposure to confounding medications. Full-thickness placental biopsies were obtained from the central region within 15 minutes of delivery, snap-frozen in liquid nitrogen, and stored at −80°C. Ethical approval was granted by the Research Ethics Committee of Ruian People’s Hospital (Approval No. YJ2024178), with written informed consent obtained from all participants.

Total RNA was isolated from placental tissues using the Tissue Total RNA Isolation Kit V2 (Vazyme Biotech, RC112-01), with purity and concentration assessed via NanoDrop spectrophotometry (Thermo Fisher Scientific; A260/A280 ratios: 1.8–2.0). Reverse transcription was performed with 1 µg RNA using HiScript III All-in-One RT SuperMix (Vazyme Biotech, R333-01) under optimized conditions: 25°C for 5 min, 50°C for 15 min, and 85°C for 5 min. qRT-PCR assays were conducted on a CFX Connect system (Bio-Rad) with Taq Pro Universal SYBR Master Mix (Vazyme Biotech, Q712-02) in 10 µL reactions (40 cycles: 95°C/10 s denaturation, 60°C/30 s annealing/extension). Melt curve analysis confirmed amplification specificity, and relative gene expression was normalized to GAPDH using the 2−ΔΔCt method. All reactions included triplicate technical replicates, with fold-change calculations relative to control samples.

All analyses were conducted within the R statistical environment (v4.4.0). Normality assumptions were verified through Shapiro-Wilk testing, with parametric comparisons (Student’s t-test) applied to normally distributed continuous variables and non-parametric alternatives (Mann-Whitney U test) for skewed distributions. Spearman’s rank correlation coefficient (ρ) quantified associations between molecular features. Unless otherwise stated, all reported p-values were two-tailed, with statistical significance defined as p < 0.05. Multiple testing correction was implemented via the Benjamini-Hochberg method for high-throughput datasets to control false discovery rate (FDR < 0.05).

Figure 1 displays the technical approach of the study, providing a concise overview of the analytical processes used in this study. The analytical flow commenced with merging transcriptomic datasets GSE75010 and GSE10588, followed by identification of DEGs. RiboRGs were intersected with DEGs to derive RiboRDEGs. Subsequent multi-modal enrichment analyses included GSEA, GO, KEGG pathway mapping, and Bayesian network inference to elucidate functional associations. WGCNA identified PE-correlated modules and key genes, while machine learning algorithms refined core diagnostic biomarkers. Model interpretability was enhanced through SHAP, with ROC curve analysis validating diagnostic efficacy using external dataset GSE10588. Immune infiltration profiling via ssGSEA revealed microenvironmental interactions of key genes. Regulatory networks encompassing mRNA-TF, mRNA-miRNA, and mRNA-RBP interactions were reconstructed to delineate molecular mechanisms. Final clinical validation confirmed differential expression patterns of candidate genes in PE cohorts.

The GSE75010 and GSE10588 datasets were merged into a combined cohort (n=200) and subjected to batch effect correction using the R package sva, followed by normalization via the R package limma. Pre-correction boxplots ( Supplementary Figure S1A ) revealed pronounced inter-batch variability in expression distributions, which resolved post-correction ( Supplementary Figure S1B ). PCA demonstrated distinct separation between original datasets along PC1 and PC2 prior to adjustment ( Supplementary Figure S1C ), whereas post-correction PCA ( Supplementary Figure S1D ) showed overlapping clusters with reduced variance contributions (PC1: 2.97%, PC2: 2.19%), confirming effective batch effect mitigation. Based on this combined cohort, we conducted differential expression analysis using the R package limma and identified 2,783 DEGs (|log2FC| > 0.1, P.adj < 0.05) between PE and control groups, comprising 1,304 upregulated and 1,479 downregulated genes ( Figure 2A ). Intersection of these DEGs with the precompiled ribosome biogenesis-related gene set (344 RiboRGs) yielded 25 RiboRDEGs, including C1QBP, DDX28, DDX51, DHX30, EXOSC2, GLUL, LSM6, MPHOSPH6, MRPL36, NCL, NOL6, PAK1IP1, POLR1B, PRKDC, RAN, RIOK1, RNASEL, RPP25, RPS27L, RRS1, SUV39H1, TBL3, TFB1M, WDR12 and XRCC5 ( Figure 2B ). The heatmap reveals the stratified clustering and expression disparity of 25 RiboRDEGs between the PE group and the control group ( Figure 2C ), which is corroborated by 3D PCA demonstrating significant intergroup separation ( Figure 2D ).

For exploring significantly dysregulated pathways in PE, we conducted GSEA analysis on the combined dataset. Results showed that top enriched terms in PE including HIF-1 signaling, AMPK signaling, and proteasome activity ( Figure 3A ). Subsequent GO and KEGG analyses of the 25 RiboRDEGs demonstrated their predominant involvement in rRNA metabolic processes (BP: ncRNA processing, ribosome biogenesis), nucleolar complexes (CC: 90S preribosome, small-subunit processome), and RNA helicase activities (MF: U3 snoRNA binding) ( Figure 3B ). KEGG pathway enrichment further implicated these genes in ribosome biogenesis and RNA degradation. Bayesian network analysis identified rRNA processing as a central hub ( Figures 3C, D ), exhibiting robust connectivity with immune response pathways (STING-mediated host immunity, IRF3-dependent IFN induction) and subcellular rRNA maturation processes.

A scale-free co-expression network was constructed using WGCNA (power = 10, scale-free fit R² = 0.9; Figure 4A ). Dynamic module detection identified 9 gene clusters ( Figures 4B, C ), with low inter-module correlations confirmed by a heatmap analysis ( Figure 4D ). The MEblue module exhibited the strongest positive correlation with PE (cor = 0.66, p = 3×10⁻²⁶; Figure 4E ). Intersection of MEblue module genes (n=955) with RiboRDEGs revealed seven hub genes (GLUL, LSM6, DDX28, NCL, RIOK1, RRS1, SUV39H1; Figure 4F ). Spearman correlation analysis demonstrated strong co-expression patterns among hub genes ( Figure 4G ), suggesting functional synergy in PE pathogenesis.

A multi-algorithm framework comprising 113 model combinations was implemented to identify optimal predictors of PE risk using seven hub genes ( Figure 5A ). Random Forest (RF) demonstrated superior diagnostic performance, achieving AUCs of 0.972 (95% CI: 0.953–0.988) in the training cohort and 0.917 (95% CI: 0.757–1.000) in the validation cohort ( Figure 5B ). The final RF model incorporated six feature genes (GLUL, DDX28, NCL, RIOK1, SUV39H1, RRS1; Table 2 ). DCA shows that this model can provide substantial clinical net benefits for clinical decision making ( Figure 5C ), while calibration curves showed strong concordance between predicted and observed PE probabilities ( Figure 5D ). Confusion matrices further validated the model’s diagnostic accuracy, achieving sensitivity >90% and specificity >85% in both training and validation cohorts ( Figure 5E ).

To elucidate the potential impact of the six feature genes on the risks associated with PE, we conducted SHAP analysis. This analysis revealed SUV39H1 as the most influential feature in the RF model (mean |SHAP| = 0.0944), followed by GLUL (|SHAP| = 0.0843), RRS1 (mean |SHAP| = 0.0649), DDX28 (mean |SHAP| = 0.0606), NCL (mean |SHAP| = 0.0534) and RIOK1 (mean |SHAP| = 0.0495) ( Figure 6A ). The expression levels of these genes have different directional effects on PE risk. Elevated SUV39H1, RRS1 and RIOK1 expression correlated with increased PE probability, whereas higher GLUL, DDX28 and NCL levels exhibited protective effects ( Figure 6B ). Interaction analysis identified synergistic risk amplification between SUV39H1 and RRS1 ( Figure 6C ), while GLUL-NCL co-expression showed concerted risk reduction ( Figure 6D ). A negative interaction between RRS1 and NCL ( Figure 6E ) suggested compensatory regulatory dynamics. Waterfall plots for representative PE (predicted probability = 0.864) and control (probability = 0.131) cases demonstrated model interpretability, with SUV39H1 contributing most substantially to risk prediction ( Figures 6F, G ).

Differential expression analysis confirmed significant downregulation of GLUL, DDX28, and NCL (p < 0.01; Figures 7A–C ) alongside upregulation of RIOK1, SUV39H1, and RRS1 (p < 0.001; Figures 7D–F ) in PE placentas. ROC analysis revealed moderate diagnostic utility for GLUL (AUC = 0.724), RIOK1 (AUC = 0.707), and SUV39H1 (AUC = 0.745) in the training cohort, with weaker performance for DDX28 (AUC = 0.613), NCL (AUC = 0.634), and RRS1 (AUC = 0.663) ( Figures 7G–L ). External validation in GSE54618 maintained moderate diagnostic accuracy for DDX28 (AUC = 0.847), RIOK1 (AUC = 0.771), SUV39H1 (AUC = 0.840), and RRS1 (AUC = 0.764), while GLUL (AUC = 0.604) and NCL (AUC = 0.688) showed limited discriminative power ( Figures 7M–R ). Protein-protein interaction (PPI) network analysis via GeneMANIA identified 20 functionally associated partners ( Figure 7S ), with co-expression (47.61%) and physical interactions (51.88%) as predominant interaction modes, suggesting that these key genes may collaboratively regulate ribosome biogenesis and RNA processing through transcriptional coordination and direct molecular binding, thereby contributing to placental dysfunction in PE pathogenesis.

To explore the systematic regulatory network of ribosome biogenesis in preeclampsia, the reconstruction of regulatory networks was conducted. Transcriptional and post-transcriptional regulatory networks were systematically mapped to elucidate molecular interactions involving key genes. The TF network, constructed using ChIPBase with a stringent filtering criterion (combined upstream/downstream supporting samples ≥8), comprised 6 key genes (GLUL, DDX28, NCL, RIOK1, SUV39H1, RRS1) interacting with 32 TFs through 59 regulatory pairs ( Supplementary Figure S2A , Supplementary Table S3 ). RBP interactions predicted via ENCORI (clipExpNum ≥10) revealed 45 functional associations between 6 key genes (GLUL, DDX28, NCL, RIOK1, SUV39H1, RRS1) and 24 RBPs ( Supplementary Figure S2B , Supplementary Table S4 ). miRNA-mediated regulation analysis identified 63 mRNA-miRNA pairs involving GLUL, SUV39H1, and NCL with 62 miRNAs (clipExpNum ≥7; Supplementary Figure S2C , Supplementary Table S5 ), highlighting transcriptional and post-transcriptional modulation of ribosome biogenesis pathways in PE.

To explore the role of Ribosome biogenesis in the immune microenvironment of PE, we conducted an ssGSEA analysis on the combined dataset. ssGSEA revealed significant dysregulation of nine immune cell types in PE placentas compared to controls (p < 0.05; Figure 8A ). Activated B cells and Th17 cells were enriched in PE, while central memory CD8+ T cells, eosinophils, immature dendritic cells, macrophages, MDSCs, memory B cells, and Th2 cells exhibited reduced infiltration. Correlation network analysis demonstrated interconnected immune cell dynamics in PE, with key genes showing cell-type-specific associations ( Figure 8B ). Notably, GLUL expression positively correlated with MDSC (r = 0.352, p < 0.001; Figure 8C ) and memory B cell abundance (r = 0.333, p < 0.001; Figure 8D ). Conversely, RRS1 displayed negative correlations with macrophages (r = -0.372, p < 0.001; Figure 8E ) and MDSCs (r = -0.337, p < 0.001; Figure 8F ), while NCL and RIOK1 inversely associated with macrophage (r = -0.334, p < 0.001; Figure 8G ) and MDSC infiltration (r = -0.301, p = 0.003; Figure 8H ), respectively.

qRT-PCR validation in placental tissues from 10 PE patients and 10 gestational age-matched controls (primers listed in Table 3 ) confirmed significant dysregulation of ribosome biogenesis-associated genes. DDX28 (p < 0.01) and RIOK1 (p < 0.001) expression was markedly elevated in PE placentas, while GLUL (p < 0.001) and NCL (p < 0.001) showed significant downregulation ( Figures 9A–D ). SUV39H1 and RRS1 exhibited non-significant expression trends ( Figures 9E, F ). Clinical assessment revealed significantly elevated systolic blood pressure in PE cases compared to gestational age-matched controls (155.1 ± 24.27 vs. 122.6 ± 7.91 mmHg, p < 0.001), with diastolic pressures similarly increased (92.9 ± 9.73 vs. 74.4 ± 3.69 mmHg, p < 0.001). No significant differences were observed in maternal age (30.4 ± 5.21 vs. 29.7 ± 4.60 years, p = 0.754), neonatal birth weight (2,724 ± 500.16 vs. 3,014 ± 503.42 g, p = 0.213), or Apgar scores at 1 min (9: 10% vs. 0%, p = 0.305) and 5 min (9: 10% vs. 0%, p = 0.305) ( Table 4 ).