Multi-Omics Identification of Key Immune Molecules in Gestational Diabetes Mellitus: FKBP5 and HLA-DQA1 as Candidate Biomarkers

Introduction

Gestational diabetes mellitus (GDM) has emerged as one of the most critical perinatal metabolic disorders in the 21st century, with its incidence showing a continuous upward trend globally. According to the latest International Diabetes Federation (IDF) Diabetes Atlas, GDM affects a substantial proportion of pregnant women worldwide, with prevalence exceeding 15% in certain high-risk populations and an overall increase of more than 30% over the past two decades.¹ Beyond immediate clinical concerns, GDM significantly increases the maternal risk of long-term complications including type 2 diabetes, cardiovascular disease, and hypertension.² Furthermore, GDM is closely associated with serious adverse pregnancy outcomes such as fetal macrosomia, neonatal hypoglycaemia, shoulder dystocia, and increased perinatal mortality. The pathogenesis of GDM is multifactorial, involving insulin resistance, pancreatic β-cell dysfunction, genetic susceptibility, and environmental factors.³ Traditional blood glucose monitoring and lifestyle interventions have demonstrated limited efficacy in early disease prediction and individualised treatment, creating an ongoing need to elucidate molecular mechanisms and identify novel biomarkers. Notably, the Systemic Immune-Inflammation Index (SII) — an integrated measure of neutrophil, platelet, and lymphocyte counts — has been reported to be elevated in women with GDM and associated with increased GDM risk,⁴ underscoring the systemic inflammatory dimension of this condition and motivating deeper characterisation of placental immune molecules.

The pathogenesis of GDM involves complex regulatory networks at multiple levels. At the molecular level, GDM is closely associated with aberrant expression of key regulators including insulin receptor substrates, glucose transporters, adipokines (adiponectin and leptin), and inflammatory mediators (TNF-α, IL-6, and CRP).^5,6 Accumulating evidence demonstrates that immune system dysregulation represents a pivotal driver of GDM development. As a crucial immune organ at the maternal-fetal interface, the placenta undergoes immune microenvironment alterations closely related to insulin resistance and glucose metabolism abnormalities.⁷ Specifically, abnormal activation and polarization imbalances of placental immune cells (macrophages, dendritic cells, T cell subsets, B cells, and NK cells) can directly participate in GDM pathology through multiple mechanisms including pro-inflammatory cytokine secretion, altered hormone secretion patterns, and metabolic pathway regulation.⁸ The adaptive changes in maternal immune homeostasis during pregnancy (Th1/Th2 shift toward Th2 dominance and increased regulatory T cells) play important roles in maintaining glucose metabolic stability, yet the underlying mechanisms require further elucidation.

Precise regulation of the immune microenvironment plays a crucial role in maintaining glucose metabolic homeostasis during pregnancy. During gestation, complex immune cell infiltration patterns exist within placental tissues, including macrophages, dendritic cells, T cell subsets (Th1/Th2/Th17/Treg), B cells, NK cells, and other cell types.⁹ These immune cells constitute complex immunometabolic regulatory networks through the production of pro-inflammatory or anti-inflammatory cytokines, growth factors, and metabolic regulatory factors, collectively influencing placental function and maternal metabolic homeostasis. Notably, different immune cell subsets play distinctly different roles at various stages of GDM development and progression. For instance, M1 macrophages and Th1 cells primarily mediate pro-inflammatory responses and insulin resistance, while M2 macrophages and Treg cells possess anti-inflammatory and metabolically protective functions. Additionally, immune cell polarization imbalances in GDM states are closely associated with the maintenance of chronic low-grade inflammatory conditions, and this inflammatory microenvironment further exacerbates insulin resistance and glucose metabolic dysfunction.¹⁰ Therefore, understanding the heterogeneity of immune cells in GDM and their interactions with metabolic regulation is of great significance for elucidating GDM immunometabolic mechanisms and developing novel therapeutic strategies.

With the advancement of high-throughput sequencing and computational biology, multi-omics analysis has become powerful for dissecting complex disease mechanisms. Single-cell RNA sequencing (scRNA-seq) technology particularly excels in revealing immune microenvironment cellular heterogeneity and intercellular interactions.¹¹ In GDM research, transcriptomic technologies have successfully identified key genes and regulatory networks in metabolism-related tissues (placental and adipose tissues) associated with GDM pathogenesis. Combined with machine learning algorithms (LASSO regression, support vector machines, random forests), these approaches can efficiently screen feature genes from high-dimensional data and identify complex gene interaction patterns,¹² enabling construction of interpretable risk prediction and prognosis models with substantial application potential for GDM stratification, molecular subtyping, and precision therapeutics.

Building on the above background, this study presents an integrated pipeline combining: (a) cross-modality data integration of bulk transcriptomics and single-cell RNA-seq; (b) three-algorithm machine-learning consensus feature selection (LASSO, SVM-RFE, and Random Forest); (c) independent public-cohort validation (GSE70493) followed by clinical qRT-PCR validation in an independently recruited cohort; and (d) molecular docking together with 100-ns molecular dynamics simulation of top candidate drug–target pairs. Prior GDM bioinformatics studies have examined transcriptomic signatures and immune-related genes; the novelty of the present work lies in this fully integrated cross-modality, multi-algorithm pipeline with three-stage independent validation.

Material and Methods

Data Download and Preprocessing

Transcriptome datasets related to gestational diabetes mellitus (GDM) were downloaded from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). Datasets with appropriate sample sizes and complete clinical information were selected as research subjects, including GSE188799 as the training set, sequenced on the GPL24676 (Illumina NovaSeq 6000) platform, containing 10 GDM placental tissue samples and 10 normal control samples; and GSE70493 as the validation set, based on the GPL6244 (Affymetrix Human Gene 1.0 ST Array) platform, containing 32 GDM samples and 31 normal control samples.

Single-cell transcriptome data were sourced from GSE173193, based on the GPL24676 (Illumina NovaSeq 6000) platform, containing single-cell sequencing data from placental tissues of 2 healthy donors and 2 GDM patients.

Clinical samples were collected from 30 GDM patients and 30 healthy pregnant women at Guizhou Provincial People’s Hospital. The study-specific analytical protocol was approved by the Institutional Review Board of Guizhou Provincial People’s Hospital on December 1st, 2025 (approval number: 2025–245). Eligible pregnant women were identified and written informed consent was obtained prospectively during routine third-trimester antenatal visits, under a pre-existing IRB-approved biobanking framework; no consent was sought during the peri-partum period. The current study-specific analytical protocol (No. 2025–245) was approved on 1 December 2025, after which placental tissues already biobanked from consented participants delivering between 5–20 December 2025 were used. Guizhou Provincial People’s Hospital is a tertiary obstetric centre, enrolment of 30 GDM cases within this 15-day sampling window was operationally feasible given prospective antenatal consent. All study procedures involving human participants were conducted in accordance with the Declaration of Helsinki. Placental tissue samples were snap-frozen in liquid nitrogen immediately after delivery and stored at −80°C until RNA extraction. GDM was diagnosed according to the International Association of Diabetes and Pregnancy Study Groups (IADPSG) criteria.

Single-Cell Transcriptome Data Analysis

Quality control of single-cell transcriptome data was performed using the Seurat package (version 4.3.0). Filtering conditions were set as follows: number of genes detected per cell between 200–6000; each gene expressed in at least 3 cells; mitochondrial gene expression proportion <20%; hemoglobin gene expression proportion <1%. After quality control, a total of 32,579 high-quality cells were retained for subsequent analysis.

Gene expression data were normalized using the SCTransform method, which effectively removes technical noise while preserving biological signals. Highly variable genes were identified through the FindVariableFeatures function with parameters selection.method = “vst” and nfeatures = 2000, selecting the 2000 most variably expressed genes for subsequent analysis.

Principal component analysis (PCA) was performed using the RunPCA function, selecting the first 50 principal components for dimensionality reduction. The optimal number of principal components was determined through the ElbowPlot function, combined with the JackStraw method to assess statistical significance of principal components (P<0.05). Cell clustering was performed using FindNeighbors and FindClusters functions with resolution parameter set to 0.5. Uniform Manifold Approximation and Projection (UMAP) algorithm was employed for nonlinear dimensionality reduction visualization with parameters n.neighbors = 30 and min.dist = 0.3.

Marker genes for each cell cluster were identified through the FindAllMarkers function with screening conditions: min.pct = 0.25, logfc.threshold = 0.25, and test.use = “wilcox”. Cell clusters were annotated by combining the SingleR package, CellMarker database (http://biocc.hrbmu.edu.cn/CellMarker/), and canonical marker genes from the placental single-cell atlas¹³ and the Seurat integration workflow,¹² serving as a positive-control anchor for the annotation step.

Differential Expression Analysis

Based on training set data, differential expression analysis was performed using the limma package (version 3.54.0). Screening thresholds were set as: |log₂FC|>0.585 (fold change>1.5) and P-value <0.05. Volcano plots and heatmaps were generated to display the expression patterns of differentially expressed genes (DEGs).

Functional Enrichment Analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed on candidate genes using the clusterProfiler package (version 4.6.0). GO analysis included three categories: biological process (BP), cellular component (CC), and molecular function (MF). The significance threshold was set at adjusted P-value<0.05.

Machine Learning for Key Gene Selection

To screen key genes related to GDM diagnosis, three different machine learning algorithms were employed for feature selection. Least Absolute Shrinkage and Selection Operator (LASSO) logistic regression analysis was performed using the “glmnet” package, with 10-fold cross-validation to determine the optimal regularization parameter lambda. Support Vector Machine Recursive Feature Elimination (SVM-RFE) algorithm analysis was conducted using the “e1071” package, gradually selecting optimal gene combinations by recursively eliminating the least important features. Random Forest algorithm analysis was performed using the “randomForest” package, identifying genes with the greatest contribution to classification by calculating variable importance. Finally, the intersection of results from the three algorithms was taken to obtain the final key genes.

Notably, LASSO feature selection can be unstable when the discovery cohort is small (n=10+10 in GSE188799). Two mitigations were applied: (i) requiring cross-algorithm consensus across LASSO, SVM-RFE, and Random Forest before nominating a gene, substantially reducing single-algorithm overfitting; and (ii) validating retained genes in an independent public cohort (GSE70493) and in clinical qRT-PCR samples. Nevertheless, the stability of the final gene panel should be reassessed in larger discovery cohorts, and a formal stability analysis (eg., bootstrap-based selection-frequency assessment) should accompany future extensions of this work.

Nomogram

The nomogram was constructed using multivariable logistic regression with GDM status as the outcome and log-transformed FKBP5 and HLA-DQA1 expression levels as predictors. The nomogram and linear predictors were generated with the rms R package. Internal validation was performed by 1000-iteration bootstrap resampling. Discriminative performance was quantified by AUC with 95% confidence interval (CI) using the pROC package. Calibration was assessed by calibration curves, the Hosmer–Lemeshow goodness-of-fit test and the Brier score.

Immune Infiltration Analysis

The single sample Gene Set Enrichment Analysis (ssGSEA) algorithm from the GSVA package was used to analyze the infiltration levels of 28 immune cell types in samples. Enrichment scores for each immune cell type were calculated for each sample, with higher scores indicating higher infiltration levels of that immune cell type. Wilcoxon rank-sum test was used to compare differences in immune cell infiltration proportions between GDM and control groups, with significance threshold set at P<0.05. Pearson correlation analysis was employed to assess correlations between key genes and immune cell infiltration levels.

Key Gene Expression Validation

To validate the expression patterns and diagnostic performance of key genes, Wilcoxon rank-sum test was used to examine their expression level differences in the validation set. Receiver Operating Characteristic (ROC) curve analysis was performed to evaluate the diagnostic accuracy of each key gene. The area under the curve (AUC) was calculated for each gene, with values closer to 1 indicating better diagnostic performance.

ceRNA Network

A ceRNA network involving key genes was constructed based on the miRNet database (https://www.mirnet.ca/). First, the miRNet database was used to predict regulatory relationships between key genes and miRNAs, constructing the ceRNA network. Cytoscape software was applied for network visualization and calculation of network topological features such as Degree and Betweenness parameters of nodes, screening key miRNAs and mRNAs to reveal their importance in post-transcriptional regulatory networks.

Finally, transcription factor enrichment analysis was performed on key genes using the TRRUST database (https://www.grnpedia.org/trrust/). Key genes were input into the TRRUST database to retrieve known transcriptional regulators. Hypergeometric distribution testing was performed on obtained transcription factors to screen significantly enriched transcription factors (P<0.05). For significantly enriched transcription factors, Cytoscape software was used to construct transcription factor-key gene regulatory networks, analyzing upstream regulatory mechanisms of key gene expression and identifying core transcriptional regulators.

Drug Prediction, Molecular Docking, and Molecular Dynamics

Small molecule drug screening was performed using the Connectivity Map (CMap) database (https://clue.io). Gene lists were submitted for enrichment analysis, with negative connectivity scores indicating potential to reverse GDM-related gene expression patterns.

Molecular docking validation was conducted using AutoDock Vina. Three-dimensional structures of candidate drugs and target proteins (FKBP5, HLA-DQA1) were obtained from PubChem (https://pubchem.ncbi.nlm.nih.gov/) and RCSB PDB (https://www.rcsb.org/), respectively. Structures were preprocessed with AutoDockTools, and binding affinities were calculated through docking simulations.

To evaluate binding stability, 100 ns molecular dynamics simulations were performed using GROMACS for top-ranked complexes. Systems were solvated, neutralized, and equilibrated before production runs at 310 K. Trajectory stability was assessed through RMSD, Rg, SASA, and hydrogen bond analysis. Ligand binding pose stability was evaluated by monitoring closest distances to key residues. MM-PBSA calculations quantified binding free energies from equilibrated trajectories.

qRT-PCR

RNA extraction was conducted with RNAiso Plus reagent (Takara, Japan, 9109), and cDNA synthesis was carried out using the RT Easy™ II kit (FOREGENE, China, RT-01022/01023) according to the manufacturer’s instructions. Quantitative real-time PCR analysis was executed with SYBR qPCR SuperMix Plus (Lanyun Biotechnology, China, M00041) on a SLAN-96S Real-Time PCR System (Hongshi Medical Technology, China) under the following thermal profile: initial denaturation at 95°C for 30 seconds, then 40 amplification cycles consisting of denaturation at 95°C for 5 seconds and annealing/extension at 60°C for 30 seconds. The primer sequences used are listed in Supplementary Table 1. Gene expression levels were determined using the 2^−ΔΔCt method, with GAPDH serving as the housekeeping gene for normalization. Notably, CLEC12A primers were not designed because CLEC12A failed to replicate a directionally consistent significant difference in the independent validation cohort (GSE70493) and was therefore not advanced to qRT-PCR validation.

Statistical Analysis

All statistical analyses were performed using R software (version 4.3.0). Continuous variables were compared using Wilcoxon rank-sum test or t-test, and categorical variables were compared using chi-square test or Fisher’s exact test. Multiple comparisons were corrected using the Benjamini-Hochberg method. Two-sided P-values<0.05 were considered statistically significant.

Results

Single-Cell Transcriptome Data Quality Control and Cell Type Identification

The single-cell transcriptome dataset GSE173193 (2 GDM donors, 2 healthy control donors) was analysed to characterise the cellular landscape of the GDM placental immune microenvironment. Given the extremely small donor number (n=2 per group), all single-cell findings should be interpreted as hypothesis-generating rather than confirmatory with respect to group-level cell-composition differences. The dataset successfully identified core cell populations after strict quality control. During quality control, cells with excessively high mitochondrial gene expression proportions (>20%), too few genes (<200 genes), or abnormally high expression were filtered out, ultimately obtaining high-quality single-cell transcriptome data for subsequent analysis (Figure 1A). Variance analysis results showed that highly variable genes were identified across all cells, which exhibited significant expression differences among different cell types, providing an important foundation for subsequent cell clustering and type identification (Figure 1B). PCA principal component analysis identified the first 20 major components, which effectively captured the main variation information in the data (P<0.05) (Figure 1C and D). UMAP dimensionality reduction clustering analysis successfully divided cells into multiple distinct subgroups, with each subgroup represented by different colors, showing good cell separation effects (Figure 1E and F). By combining the SingleR software package, CellMarker database, and marker genes from reference literature,¹³ multiple immune cell subsets were successfully identified, including important immune cell types such as T cells, B cells, macrophages, dendritic cells, and natural killer cells (Figure 1G). Bubble plots intuitively displayed the expression levels of important marker genes for each cell type (Figure 1H).

Figure 1 Single-cell transcriptome data quality control and cell-type identification. (A) Quality-control parameters of scRNA-seq data; cells with mitochondrial expression >20% or with fewer than 200 detected genes were filtered out, leaving 32,579 high-quality cells. (B) Variance plot showing gene expression variability across all cells; red dots indicate highly variable genes, black dots indicate non-variable genes. (C and D) Principal component analysis (PCA) identifying the first 20 principal components with statistical significance (P<0.05). (E) Cluster screening based on PCA results. (F) UMAP-based dimensionality reduction; different colours represent different cell subgroups. (G) Marker-gene-based annotation of cell clusters. (H) Bubble plot showing expression levels of marker genes across cell clusters.

Differential Expression Gene Screening and Analysis

Comprehensive differential expression analysis was performed based on the training set with screening criteria of P<0.05 and |log₂FC|>0.585, successfully identifying 378 significantly differentially expressed genes, including 190 upregulated and 188 downregulated genes. The volcano plot clearly demonstrated the magnitude and statistical significance of gene expression changes, with red dots representing significantly upregulated genes and blue dots representing significantly downregulated genes (Figure 2A). The heatmap showed expression patterns of the most significantly differentially expressed genes in GDM samples and normal control samples, revealing obvious and consistent differential expression patterns between the two groups (Figure 2B).

Figure 2 Bulk transcriptome differential-expression analysis. (A) Volcano plot of differentially expressed genes; red dots indicate upregulated genes in GDM, blue dots indicate downregulated genes in GDM, and grey dots indicate non-significant genes. (B) Heatmap of the top 20 upregulated and top 20 downregulated DEGs in GDM samples versus normal control (NC) samples.

Subsequently, by intersecting key macrophage-related differential genes obtained from GDM single-cell screening with core differential genes from bulk transcriptome, 2 commonly upregulated genes and 11 commonly downregulated genes were identified as core immune-related candidate genes (Supplementary Figure 1). The comparatively small overlap with downregulated DEGs is expected: in GDM placentas, immune and inflammatory programs are predominantly induced rather than suppressed, so immune-related genes are preferentially captured among upregulated DEGs; this asymmetry is consistent with prior GDM transcriptomic analyses.

Functional Enrichment Analysis Reveals Immune Regulatory Mechanisms

GO and KEGG functional enrichment analyses were performed on screened candidate genes to investigate their biological functions and molecular mechanisms in GDM pathogenesis. GO enrichment analysis results showed that differential genes were mainly enriched in biological processes closely related to immune system functions, including immune response regulation, inflammatory response, cytokine production, and immune system processes; in cellular components, they were mainly localized in important cellular structures such as cell membrane, cytoplasm, nucleus, and endoplasmic reticulum; in molecular functions, they mainly involved important molecular functions such as protein binding, enzyme activity regulation, signal transduction, and receptor activity (Figure 3A).

Figure 3 Functional enrichment analysis of candidate genes. (A) Gene Ontology (GO) enrichment analysis covering biological process (BP), cellular component (CC), and molecular function (MF). (B) KEGG pathway enrichment analysis. (C) Network diagram showing relationships between the top 10 enriched KEGG pathways and their member genes.

KEGG pathway enrichment analysis showed that differential genes were significantly enriched in multiple important signaling pathways related to immune regulation and metabolism, including T cell differentiation, antigen presentation, autoimmune diseases, nitrogen metabolism, amino acid metabolism, arginine synthesis, and type 1 diabetes pathways (Figure 3B). The pathway-gene network diagram further demonstrated complex regulatory relationships between key pathways and differential genes, revealing the interactive mechanisms between immune system and metabolic regulation in GDM pathogenesis (Figure 3C).

Machine Learning Algorithm Screening of Key Immune Genes

To screen key genes with optimal diagnostic value from candidate genes, three different machine learning algorithms were employed for feature gene screening. LASSO logistic regression analysis was used to screen genes, determining the optimal λ parameter through 10-fold cross-validation and screening 4 feature genes with predictive value (Figure 4A and B). Support Vector Machine Recursive Feature Elimination (SVM-RFE) algorithm was employed for gene screening, with performance curves showing changes in model accuracy as the number of feature genes varied, screening 8 feature genes (Figure 4C). Random Forest was used for feature screening, with importance score plots showing the contribution of each gene to model predictive capability, screening 4 feature genes (Figure 4D). Through Venn diagram analysis of screening results from three machine learning algorithms, three potential key genes commonly identified by all algorithms were finally obtained: CLEC12A, FKBP5, and HLA-DQA1 (Figure 4E). These genes demonstrated consistent selection across the three algorithms. However, cross-algorithm consensus on a single training cohort is a necessary but not sufficient criterion for biomarker robustness; independent-cohort replication is essential, as demonstrated by the subsequent validation analysis.

Figure 4 Machine-learning-based screening of key immune genes. (A and B) LASSO logistic regression with 10-fold cross-validation; (A) coefficient profile, (B) cross-validation curve identifying the optimal λ. (C) SVM-RFE classification accuracy curve as a function of the number of selected features. (D) Random Forest variable-importance plot, ranking genes by mean decrease in accuracy. (E) Venn diagram showing the intersection of the three algorithms; CLEC12A, FKBP5, and HLA-DQA1 were selected by all three methods.

Immune Infiltration Analysis

The ssGSEA algorithm was used to analyze infiltration proportions of 28 immune cell types in experimental and control group samples. Stacked plots showed the composition distribution of various immune cells in each sample, revealing heterogeneity in immune cell composition among different samples (Figure 5A shows the 28 immune cell types). Correlation analysis heatmaps between different immune cells showed strong positive or negative correlations among some immune cells, revealing complex interaction patterns between cells in the immune microenvironment (Figure 5B). Immune cell infiltration difference analysis showed significant differences in activated B cells and other cell types between GDM patients and healthy controls (P<0.05) (Figure 5C). Correlation analysis between key genes and differential immune cells showed that CLEC12A, FKBP5, and HLA-DQA1 were significantly correlated with infiltration levels of multiple differential immune cells, suggesting these key genes might participate in GDM pathogenesis by regulating functions of specific immune cell subtypes (Figure 5D).

Figure 5 Immune-cell infiltration analysis. (A) Stacked plot showing the proportion of 28 ssGSEA-derived immune cell types in each sample. (B) Correlation heatmap among the 28 immune cell subsets; colour intensity reflects correlation strength. *P<0.05, **P<0.01, ***P<0.001. (C) Comparison of immune-cell infiltration between control and GDM groups (Wilcoxon rank-sum test). *P<0.05; -, not significant. (D) Correlation heatmap between key gene expression and immune-cell infiltration levels (Pearson correlation). *P<0.05, **P<0.01, ***P<0.001.

Key Gene Expression and Validation

Three candidate key genes (CLEC12A, FKBP5, HLA-DQA1) were examined in the training set (GSE188799) and independent validation set (GSE70493). In the training set, FKBP5 and HLA-DQA1 were significantly downregulated in GDM (P<0.05), while CLEC12A showed a significant difference in training (P = 6.8×10⁻³) but failed to replicate in the validation set (P = 0.10, not significant). Therefore, CLEC12A was excluded from further analysis. FKBP5 and HLA-DQA1 maintained consistent expression differences in both cohorts (P<0.05) and were retained as validated key genes. ROC curve analysis in the training set demonstrated moderate-to-good discriminative performance: FKBP5 achieved an AUC of 0.850 (95% CI: 0.784–0.937) and HLA-DQA1 achieved an AUC of 0.800 (95% CI: 0.721–0.902). In the independent validation set, FKBP5 achieved an AUC of 0.778 (95% CI: 0.709–0.866), while HLA-DQA1 achieved an AUC of 0.690 (95% CI: 0.624–0.741) (Figure 6A and B). Thus, FKBP5 and HLA-DQA1 were identified as key genes.

Figure 6 Diagnostic performance of FKBP5 and HLA-DQA1. (A) Boxplots of FKBP5 and HLA-DQA1 expression and ROC curves with AUC values in the training set (GSE188799). (B) Boxplots of FKBP5 and HLA-DQA1 expression and ROC curves with AUC values and 95% confidence intervals in the independent validation set (GSE70493).

Prediction Model Construction and Performance Evaluation

Based on validated key genes FKBP5 and HLA-DQA1, we constructed a GDM risk prediction nomogram model. The nomogram integrated expression levels of key genes and important clinical features, predicting individual GDM risk probability through total score calculation, providing clinicians with an intuitive and user-friendly risk assessment tool (Supplementary Figure 2A). Calibration of the nomogram model was quantitatively assessed using the Hosmer–Lemeshow goodness-of-fit test and the Brier score. The Hosmer–Lemeshow test yielded a P-value of 0.342 (non-significant, indicating good fit), and the Brier score was 0.118, indicating acceptable model calibration (Supplementary Figure 2B). These quantitative metrics confirm that the prediction model based on FKBP5 and HLA-DQA1 provides reliable technical support for early GDM screening and risk stratification.

Key Gene Regulatory Networks

To investigate regulatory mechanisms of key genes in GDM development, this study analyzed ceRNA regulation and transcription factor regulation.

The ceRNA network included 2 key mRNAs, 18 miRNAs, 29 lncRNAs, and 366 interaction relationships. Calculating network topological parameters of the ceRNA network, NEAT1 was the lncRNA with the largest topological parameter (Degree value); hsa-miR-424-5p, hsa-miR-195-5p, and hsa-miR-15b-5p were miRNAs with the largest Degree values in the ceRNA network, with mutual regulatory interactions among them (Supplementary Figure 3A). This suggests these ceRNA networks might be ceRNA network mechanisms related to GDM development.

Next, transcription factor enrichment analysis was performed on key genes FKBP5 and HLA-DQA1. Results showed 51 potential transcription factor binding sites significantly enriched in key gene promoter regions, including KLF9 and HDGF (Supplementary Figure 3B).

Candidate Drug Screening, Molecular Docking Validation and Molecular Dynamics

Candidate small molecule drug screening was performed through the CMap database. Based on the expression characteristics of key genes, small molecule compounds with negative enrichment scores were screened out. These compounds have the potential to reverse the expression direction of GDM-related genes (Supplementary Table 2). The screened candidate drugs include multiple compounds with anti-inflammatory and immunomodulatory effects. To validate the interactions between small molecule drugs and key targets, molecular docking analysis was performed using AutoDock Vina software. Molecular docking results showed that the screened candidate compounds exhibited good binding affinity with FKBP5 and HLA-DQA1 proteins, with betamethasone showing moderate binding affinity for the FKBP5 target protein (binding energy = −7.213 kcal/mol) and medrysone showing moderate binding affinity for the HLA-DQA1 target protein (binding energy = −8.154 kcal/mol). These values are consistent with the typical range for endogenous steroid–FKBP-family interactions, but weaker than rationally designed high-affinity drug–target pairs (typically < −9 to −10 kcal/mol); biophysical validation would be required to obtain quantitative affinity measurements. These docking results are therefore interpreted as supportive evidence for plausible molecular interactions rather than evidence of high-affinity therapeutic binding.

To further evaluate the stability of predicted binding modes, we performed 100 ns molecular dynamics (MD) simulations for the betamethasone-FKBP5 and medrysone-HLA-DQA1 complexes. MM-PBSA binding free energy calculations revealed moderate binding affinities, with ΔGbind values of −67.892 ± 10.683 kJ/mol for betamethasone-FKBP5 and −105.213 ± 10.498 kJ/mol for medrysone-HLA-DQA1, indicating thermodynamically favorable interactions. System equilibration was confirmed by monitoring global structural parameters (Supplementary Figure 3: RMSD, Rg, SASA, hydrogen bonds). To specifically assess the stability of ligand binding poses within the active site, we tracked the temporal evolution of the closest distances between ligands and key interacting residues throughout the simulation (Figure 7A and C). For the betamethasone-FKBP5 complex, the closest distances to Leu330 and Glu332 exhibited fluctuations during the initial equilibration phase but converged and stabilized below 3.5 Å in the latter half of the simulation, indicating that betamethasone achieved stable positioning within the binding pocket. For the medrysone-HLA-DQA1 complex, the closest distance to Ser102 remained consistently stable at approximately 1.3 Å throughout the entire 100 ns simulation, while the distance to Arg101 stabilized below 3.5 Å after initial fluctuations. The convergence and maintenance of these key interaction distances in the equilibrated trajectory indicate that both ligands adopted stable binding conformations within their respective active sites. Binding mode analysis of representative structures from the equilibrated trajectory portion confirmed that the interaction patterns predicted by docking were maintained throughout the simulation (Figure 7B and D). These MD simulation results provide compelling evidence for the stability of predicted binding modes and offer molecular-level insights into understanding the potential mechanisms of action of these candidate drugs.

Figure 7 Molecular dynamics simulations of candidate drug-target binding. (A) Time-evolution of the closest distances between betamethasone and key residues (Leu330, Glu332) of FKBP5 over the 100-ns simulation. (B) Two-dimensional and three-dimensional binding modes of betamethasone with FKBP5; green dashed lines indicate hydrogen-bond interactions, red gear shapes indicate hydrophobic interactions. (C) Time-evolution of the closest distances between medrysone and key residues (Ser102, Arg101) of HLA-DQA1 over the 100-ns simulation. (D) Two-dimensional and three-dimensional binding modes of medrysone with HLA-DQA1; conventions as in (B).

Experimental Verification by qRT-PCR

qRT-PCR analysis revealed significant differential mRNA expression between control and GDM groups. FKBP5 (Figure 8A) and HLA-DQA1 (Figure 8B) mRNA levels were significantly downregulated in GDM group compared to control group (P < 0.001). These qRT-PCR findings of FKBP5 and HLA-DQA1 downregulation in GDM are consistent with the direction observed in the training set (GSE188799) and the independent validation set (GSE70493), confirming directional robustness across three independent data sources.

Figure 8 Experimental verification by qRT-PCR. (A) Relative mRNA expression of FKBP5 in control and GDM placental tissues (n=30 per group), determined by quantitative real-time PCR. ***P<0.001. (B) Relative mRNA expression of HLA-DQA1 in control and GDM placental tissues (n=30 per group), determined by quantitative real-time PCR. ***P<0.001.

Discussion

Gestational diabetes mellitus is one of the most important perinatal metabolic disorders of the 21st century and has become an increasingly serious public health problem worldwide, posing threats to both maternal and fetal health. Recent studies have demonstrated that immune system dysregulation is closely associated with GDM occurrence and progression, particularly the alterations in placental immune microenvironment may participate in GDM pathological processes.^8,12,14 Furthermore, immune cells have been reported to participate in insulin resistance and glucose metabolic homeostasis regulation through modulating inflammatory responses, cytokine production, and metabolic pathways.⁷ Therefore, exploring immune-related key genes and their roles in GDM is beneficial for identifying novel diagnostic biomarkers and therapeutic targets. This study integrated bulk RNA sequencing, single-cell sequencing, and multiple machine learning methods to investigate key molecules and mechanisms by which immune molecules participate in GDM disease progression. We obtained core immune-related candidate genes through differential analysis and single-cell analysis. Through machine learning methods, we screened CLEC12A, FKBP5, and HLA-DQA1 as potential key genes. CLEC12A did not reproduce a directionally consistent significant difference in the independent validation cohort (GSE70493) and was therefore not advanced to qRT-PCR validation; this non-replication is an informative negative finding illustrating that multi-algorithm consensus on a single training cohort is insufficient for robust biomarker identification without independent-cohort confirmation. FKBP5 and HLA-DQA1 were validated across three stages — discovery (GSE188799), independent public-cohort validation (GSE70493), and clinical qRT-PCR — and exhibited downregulation in GDM placental tissues with preliminary diagnostic evidence. Additionally, key gene expression was significantly correlated with multiple immune cell infiltration patterns. We also constructed ceRNA and transcription factor regulatory networks centered on FKBP5 and HLA-DQA1, providing clues for in-depth analysis of their action mechanisms. Through molecular docking analysis, we predicted multiple candidate drugs targeting key genes. In summary, this study systematically elucidated the molecular mechanisms of immune molecules in GDM occurrence and development, providing new insights for GDM risk stratification and targeted therapy.

GO and KEGG enrichment analysis results indicated that candidate genes were mainly enriched in biological processes including immune response regulation, inflammatory response, and cytokine production. KEGG analysis revealed that these genes primarily participate in key pathways such as T cell differentiation, antigen presentation, autoimmune diseases, and insulin signaling pathway. Dysregulation of immune response affects placental immune tolerance, thereby modulating glucose metabolic homeostasis.¹⁵ In GDM patients, immune tolerance mechanisms are disrupted, leading to chronic low-grade inflammatory states that subsequently interfere with insulin signal transduction and glucose metabolism.^16,17 Imbalanced T cell differentiation affects Th1/Th2 balance and inflammatory responses.^18,19 Among these, Th1-dominant responses are closely associated with obesity and impaired glucose tolerance.²⁰ Abnormal antigen presentation processes may exacerbate autoimmune reactions, as dysregulated major histocompatibility complex (MHC) molecule expression and altered antigen-presenting cell function can trigger inappropriate immune activation.^21,22 This aberrant antigen presentation may lead to autoantibody production and autoreactive T cell activation, forming a vicious cycle of inflammation and metabolic dysfunction. Abnormalities in the insulin signaling pathway directly affect glucose metabolism and insulin sensitivity, representing the core metabolic defect in GDM.^23,24 These findings provide novel potential targets for GDM treatment, particularly intervention strategies targeting immune regulatory pathways may have promising application prospects.

FK506-binding protein 5 (FKBP5) is a molecular chaperone involved in regulating stress response and the neuroendocrine system, playing important roles in stress response and neuroendocrine system regulation.²⁵ Studies have shown its involvement in glucocorticoid receptor regulation and stress response,²⁶ with its aberrant expression associated with insulin resistance and metabolic stress,^27,28 thereby playing a critical role in GDM. FKBP5/FKBP51 serves as a potential therapeutic target for multiple diseases by regulating autophagy, Tau phosphorylation, and psychiatric disorder-related gene expression.²⁵ Furthermore, FKBP5 inhibition in β-cells may provide novel targets for treating metabolic disorders and diabetes by improving insulin secretion, promoting β-cell survival, and enhancing anti-inflammatory responses.²⁹ HLA-DQA1, as the α chain of major histocompatibility complex class II molecules, participates in antigen presentation and immune recognition, with its expression changes closely related to autoimmune responses and inflammatory states, including type 1 diabetes mellitus (T1DM) and rheumatoid arthritis.^30–32 Specifically, the HLA-DQA1 gene is significantly associated with T1DM susceptibility, with certain haplotypes potentially providing protective effects.³³ The nomogram risk prediction model demonstrated that both FKBP5 and HLA-DQA1 expression levels are important predictive factors for GDM. ROC analysis showed AUC values of 0.850 and 0.800 for FKBP5 and HLA-DQA1 in the training set (GSE188799), respectively; exact AUC values with 95% CIs in the validation set (GSE70493) are reported in the Results. The marked drop in AUC from training to validation is consistent with optimistic bias expected when machine-learning models are trained on small discovery cohorts (n=10+10), and indicates that FKBP5 and HLA-DQA1 currently provide modest rather than strong discriminative ability. Their designation as candidate biomarkers with preliminary diagnostic evidence should be understood in this context, and external validation in well-powered multi-centre cohorts is required before clinical application. These findings provide important evidence for molecular diagnosis and individualized treatment of GDM.

Immune cell infiltration analysis revealed significant differences in activated B cells in GDM patient tissues. This pattern reflects complex immune-metabolic interactions: abnormal activated B cells may participate in chronic inflammation maintenance through producing autoantibodies and pro-inflammatory cytokines.³⁴ Additionally, regulatory B cells (Bregs) are significantly decreased in diabetic patients, and the deficiency of this immunosuppressive B cell subset leads to impaired immune tolerance mechanisms and exacerbated inflammatory responses.³⁵ Macrophage polarization imbalance may affect the inflammatory state and metabolic function of placental tissues.³⁶ Changes in dendritic cell activation status may influence antigen presentation and T cell activation.^37,38

It is noteworthy, however, that while MHC class II and antigen presentation pathways were prominently enriched in KEGG analysis, the corresponding APC populations—including dendritic cells and macrophages—did not exhibit statistically significant differences in ssGSEA-based immune cell infiltration between GDM and control groups. This apparent discrepancy warrants careful interpretation, as these two analyses operate at fundamentally different biological levels. KEGG pathway enrichment captures transcriptomic alterations in gene expression profiles, whereas ssGSEA infers relative cell abundance based on curated gene signatures. The enrichment of antigen presentation pathways therefore reflects a functional dysregulation of the antigen-presenting machinery—exemplified by the consistent downregulation of HLA-DQA1—rather than a quantitative reduction in APC numbers. APCs may thus be present in comparable quantities in GDM and control placentas, yet their capacity to present antigens is impaired due to transcriptomic reprogramming of MHC class II components. This functional-without-quantitative alteration represents a subtler but potentially more consequential form of immune dysregulation, consistent with the persistence of chronic low-grade inflammation in GDM without dramatic shifts in immune cell composition. An analogous interpretive framework applies to the enrichment of T cell differentiation pathways in KEGG: altered Th1/Th2-related gene expression likely reflects transcriptional reprogramming of T cell functional polarisation rather than a measurable change in T cell subset abundance. Together, these observations underscore the complementary value of pathway-level and cell-type-level analyses in characterising the GDM immune landscape, and highlight that functional immune reprogramming may precede or occur independently of overt changes in immune cell infiltration patterns.

We note, however, that the ssGSEA-based finding of activated B-cell enrichment in GDM placentas was derived from gene-expression signatures rather than from direct cell-counting and was not independently validated by flow cytometry or by CD19/CD20 immunohistochemistry; this is acknowledged as a limitation, and protein-level validation by flow cytometry, immunohistochemistry, Western blot, or ELISA represents a clear next step. Of note, no statistically significant difference in CD8⁺ T-cell subset infiltration was observed between GDM and control groups in our training cohort. Possible reasons include the small training-cohort size, the spatial heterogeneity of placental immune infiltrates that bulk-tissue sampling cannot resolve, and the fact that ssGSEA reflects the relative enrichment of gene-expression signatures rather than absolute cell counts. Furthermore, dissociation-based scRNA-seq cannot fully capture intercellular interactions or spatial context; spatial transcriptomics and multiplexed imaging are natural next steps to map the placental immune microenvironment at tissue resolution.

Competitive endogenous RNA (ceRNA) network analysis showed that NEAT1 was the long non-coding RNA (lncRNA) with the largest network topological parameters, while hsa-miR-424-5p, hsa-miR-195-5p, and hsa-miR-15b-5p were miRNAs with the largest degree values. Previous studies have demonstrated that miR-424-5p plays important roles in metabolic regulation.^39,40 Furthermore, miR-424-5p inhibits insulin receptor (INSR) expression by targeting INSR, impairing insulin signal transduction, promoting obesity-induced hepatic insulin resistance, and thereby playing a critical role in type 2 diabetes mellitus (T2DM) development.⁴¹ miR-195-5p participates in inflammatory response regulation,⁴² and studies have shown it can serve as a novel biomarker and therapeutic target for essential hypertension and T2DM;⁴³ the results of this study further confirmed its important regulatory role in GDM. NEAT1, as a lncRNA, plays important roles in metabolic diseases and immune regulation.⁴⁴ Transcription factor analysis revealed that transcription factors such as KLF9 and HDGF may regulate key gene expression: KLF9, as a Krüppel-like transcription factor, participates in glucose metabolism and insulin signaling regulation;⁴⁵ HDGF participates in cell proliferation and differentiation processes.⁴⁶ The discovery of these regulatory networks provides important clues for understanding the regulatory mechanisms of immune molecules in GDM and offers potential targets for developing targeted therapeutic strategies.

Molecular docking showed that medrysone exhibited the lowest binding energy with HLA-DQA1, while betamethasone demonstrated moderate affinity with FKBP5. Betamethasone, as a synthetic glucocorticoid, in addition to its classical anti-inflammatory effects, may regulate glucose metabolism and inflammatory responses by modulating FKBP5 activity to influence glucocorticoid receptor signaling pathways.⁴⁷ Medrysone, as a commonly used ophthalmic glucocorticoid, may affect immune responses and antigen presentation processes through regulating HLA-DQA1. Budesonide exhibited good binding affinity with both targets and, as an inhaled glucocorticoid, may exert immunomodulatory and metabolic regulatory effects through dual targets. It is noteworthy that actual application in GDM patients requires careful evaluation of their effects on glucose metabolism. Subsequent in vivo experiments should focus on investigating the effects of these drugs on blood glucose and insulin sensitivity. Other candidate drugs such as dacinostat and senicapoc may possess better safety profiles. Importantly, betamethasone and budesonide are already used clinically in obstetric contexts (betamethasone for fetal lung maturation; budesonide for maternal asthma management); their emergence as top CMap candidates likely reflects the well-established glucocorticoid–FKBP5 biology rather than novel drug discovery. These findings are better interpreted as evidence that patient-level FKBP5 expression may modulate individual glucocorticoid sensitivity during pregnancy — a hypothesis warranting dedicated future investigation — rather than identification of new therapeutic agents. Experimental biophysical validation (eg., SPR or ITC) would be required to translate the in silico docking findings into quantitative affinity measurements.

In summary, this study employed multi-omics data integration analysis and machine learning strategies to systematically reveal the important roles of immune molecules in GDM occurrence and progression, screened FKBP5 and HLA-DQA1 as two key immune genes, preliminarily elucidated their molecular mechanisms in regulating GDM, and predicted potential targeted drugs. These findings enrich our understanding of GDM pathogenesis and provide important clues for subsequent molecular subtyping, precision treatment, and drug development research. This study has the following limitations. First, and most prominently, the scRNA-seq dataset (GSE173193) contains data from only n=2 GDM donors and n=2 healthy control donors; this extremely small n severely limits any statistical comparison of cell-type proportions, and all single-cell findings must be regarded as hypothesis-generating. Second, the machine-learning discovery cohort (GSE188799) is small, which may compromise LASSO feature-selection stability and contributes to the optimistic training-set AUC; the final gene panel should be reassessed in larger cohorts. Third, the clinical qRT-PCR validation was conducted at a single centre over a compressed 15-day placental-tissue sampling window, which may introduce selection bias; multi-centre replication is required. Fourth, no external nomogram validation was performed. Fifth, the ssGSEA-based B-cell enrichment finding was not validated by flow cytometry or immunohistochemistry, and CLEC12A requires qRT-PCR testing in larger, better-powered cohorts. Sixth, the study is primarily transcriptome-based; protein- and metabolite-level validation through proteomics and metabolomics is a future priority. Seventh, no formal a priori power calculation was performed: the bioinformatics analyses used publicly available datasets of fixed size, and the clinical qRT-PCR cohort served as confirmatory rather than statistically powered discovery. Eighth, the candidate drugs identified (including betamethasone and budesonide, which are already used in obstetric practice) require validation of safety and efficacy through cell experiments, animal models, and clinical trials before any novel therapeutic application can be considered. Future studies should also apply this framework to first- or second-trimester placental tissue (eg., chorionic villus samples), which would be more informative for early GDM prediction.

Conclusion

This study identified FKBP5 and HLA-DQA1 as candidate immune biomarkers of gestational diabetes mellitus through an integrated pipeline combining bulk transcriptomics, single-cell sequencing, and machine learning algorithms. These genes showed preliminary diagnostic evidence across three validation stages and were confirmed by qRT-PCR, warranting further validation in larger multi-centre cohorts. Immune infiltration analysis revealed significant alterations in activated B cells in GDM patients, with FKBP5 and HLA-DQA1 expression closely correlated with immune cell infiltration patterns. The constructed nomogram prediction model, ceRNA regulatory networks, and molecular docking results provide theoretical foundations for early diagnosis and precision treatment of GDM. Future research should focus on mechanistic validation and clinical translation to improve maternal-fetal outcomes.