Four transcriptomic datasets—GSE48424, GSE75010, GSE39653, and GSE98793—were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). The GSE75010 dataset, comprising placental transcriptome data related to preeclampsia (PE), includes 77 control samples and 80 placental tissue samples from PE patients. The peripheral blood transcriptome dataset GSE48424 includes 18 control samples and 18 samples from PE patients. Additionally, GSE39653, containing PBMC transcriptomic data from patients with DD, was employed as the training set, consisting of 24 control samples and 29 samples from DD patients. The GSE98793 dataset, derived from microarray data of peripheral blood from DD patients, was used for validation. To ensure data comparability, gene expression matrices across all datasets were normalized via the “normalizeBetweenArrays” function in R. Differentially expressed genes (DEGs) were subsequently identified via the limma package, with selection thresholds set at FDR < 0.05 and |logFC| > 0.233. The resulting DEGs were visualized via the pheatmap and ggplot2 packages to generate heatmaps and volcano plots that display the gene expression characteristics. A total of 3,946 human genes encoding secreted proteins were retrieved from the Human Protein Atlas database (https://www.proteinatlas.org/).

To identify gene modules strongly associated with the phenotype and to uncover hub genes within the network, weighted gene coexpression network analysis (WGCNA) was conducted. On the basis of the peripheral blood transcriptomic dataset GSE39653 from individuals with DD, a gene coexpression network was constructed via the “WGCNA” package in R. The 25% most highly variable genes among the differentially expressed genes were selected as inputs for network construction. A soft-thresholding power of β = 4 and a scale-free topology fit index of R² = 0.9 were applied to ensure that the resulting network conformed to a scale-free topology. Hierarchical clustering was then performed via the topological overlap matrix (TOM) to generate a gene dendrogram, with individual modules distinguished by unique colors. The most significantly associated module was selected for downstream analysis.

To investigate protein–protein interaction (PPI) relationships between PE-associated secreted proteins and key genes implicated in DD, a PPI network was constructed via the STRING database (version 12.0, https://www.string-db.org), with an interaction confidence score threshold set at 0.4 (28). The resulting network was visualized via Cytoscape software (version 3.8.2), and functional module identification was performed via the MCODE plugin. The MCODE parameters were configured as follows: degree cut-off = 2, node score cut-off = 0.2, k-core = 2, and maximum depth = 100. The three functional modules with the highest scores were selected for downstream bioinformatic analyses.

To elucidate the functional characteristics and biological pathways associated with candidate genes related to PE and DD, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed via the “clusterProfiler” and “enrichplot” packages in R. Gene annotation data were obtained from the “org.Hs.eg.db” file available in the Bioconductor database. The enrichment results were visualized via the “ggplot2” package.

To identify potential diagnostic biomarkers for DD associated with PE, two machine learning approaches were employed for gene selection. Initially, feature selection was performed via least absolute shrinkage and selection operator (LASSO) logistic regression implemented via the “glmnet” package in R (29). The gene importance was subsequently evaluated via the random forest (RF) algorithm via the “randomForest” package (30). Genes identified by both LASSO regression and RF analysis were intersected and considered candidate diagnostic markers for PE-associated DD.

To assess the diagnostic efficacy of candidate genes for PE-related DD, a multivariate logistic regression model was established by integrating key hub genes via the “Irm” package in R to predict the risk of disease. A nomogram was generated to intuitively display the model and variable contributions (31). The model’s discrimination ability was evaluated via receiver operating characteristic (ROC) curves, and the area under the curve (AUC) was calculated. Calibration curves were used to assess the consistency between the predicted and observed probabilities. Decision curve analysis (DCA) was further applied to quantify the clinical benefit of the model. External validation of the nomogram’s predictive accuracy was performed using the independent DD cohort GSE98793.

Immune cells exhibit distinct patterns of infiltration and residency during disease initiation and progression, offering crucial insights into their roles in pathophysiological mechanisms. On the basis of the GSE39653 dataset of DD, the CIBERSORT algorithm was employed to quantify the relative infiltration levels of immune cells within disease tissues. The Wilcoxon rank-sum test was used to assess the differences in the proportions of 22 immune cell types between the DD group and healthy controls.FDR correction was applied to adjust for multiple comparisons. Spearman’s rank correlation analysis was conducted to evaluate the relationship between the expression levels of diagnostic biomarkers and immune cell infiltration. Furthermore, gene set enrichment analysis (GSEA) of characteristic genes was performed via the “org.Hs.eg.db” and “clusterProfiler” packages. Significantly enriched signalling pathways in DD were visualized via the “enrichplot” package.

Putative transcription factors (TFs) associated with the hub genes were predicted via the KnockTF v2.0 database (https://bio.liclab.net/KnockTFv2/index.php). To visualize the regulatory mechanisms of the identified hub genes, a TF–hub gene regulatory network was constructed via Cytoscape software.

To investigate potential environmental exposures associated with the characteristic genes, chemical compounds linked to these genes were retrieved from the Comparative Toxicogenomics Database (CTD, https://ctdbase.org/), with a particular focus on phthalates. A gene–chemical interaction network was constructed via Cytoscape to delineate potential roles of the identified genes within toxicological networks. To further explore the molecular interaction mechanisms between phthalates and key hub gene products, molecular docking simulations were performed via AutoDock Vina, and three-dimensional visualization of docking conformations was conducted via PyMOL. Protein structural data for target genes were obtained from the Protein Data Bank (PDB), while ligand structures were retrieved from the PubChem database.

R software (version 4.4.1) and various packages were used for statistical analysis and visualization. A p value of less than 0.05 was considered statistically significant.

The analytical workflow of this study is illustrated in Figure 1. Two raw datasets, comprising samples from PE patients and healthy controls, were downloaded from the GEO database. After normalization of the raw microarray data, differential expression analysis was performed via the “limma” package. The results revealed a total of 1,626 differentially expressed genes (DEGs) in placental tissue, including 781 upregulated and 845 downregulated genes. In peripheral blood samples, 270 DEGs were detected, comprising 130 upregulated and 140 downregulated genes. The expression patterns of DEGs in placental tissue and peripheral blood from PE patients are visually presented via volcano plots and heatmaps (Figures 2A, B). Furthermore, to investigate PE-associated secretory proteins, a set of 3,946 human genes encoding secreted proteins was intersected with the PE DEGs. This analysis revealed 434 secretory protein genes related to PE, including 400 derived from placental tissue and 37 from peripheral blood samples (Figure 2C).

Weighted gene coexpression network analysis (WGCNA) was performed on the depression dataset GSE39653 to identify the gene modules most closely associated with depression. A soft-thresholding power of β=10 was chosen to ensure that the constructed coexpression network conformed to scale-free topology and met the average connectivity criteria (Figure 3A). The minimum module size was set to 150, and a module merging threshold of 0.25 was applied. Hierarchical clustering dendrograms for depression-related modules were generated (Figure 3B). Seven coexpression modules, each assigned a distinct color, were identified. Correlation analysis revealed that the blue module presented the strongest positive correlation with DD (r = 0.23, p = 0.01) (Figure 3C). Additionally, 19 gene modules were identified within the depression dataset, among which the MEblack module showed the strongest negative correlation (Figure 3D). Consequently, 1,165 key genes from the gray module were selected for subsequent analysis (Figure 3D).

Clinical studies have revealed that patients with preeclampsia (PE) have a greater risk of developing depressive disorder (DD), suggesting a potential causal link between these two conditions. To investigate this association, a protein–protein interaction (PPI) network comprising 434 PE-associated secreted proteins and 1,165 MDD-related pathogenic genes was constructed via the STRING database. The PPI network was further analysed via the MCODE plugin in Cytoscape, which identified 24 significant modules. Among these, the top three modules based on MCODE scores contained a total of 105 genes, which were considered candidate DD pathogenic genes associated with PE (Figure 4A). GO enrichment analysis indicated that these genes were predominantly involved in biological processes such as leukocyte migration, cell chemotaxis, and positive regulation of cytokine production. These genes were localized to cellular components, including the collagen-containing extracellular matrix and the endoplasmic reticulum lumen, and were functionally associated with cytokine activity and chemokine activity (Figure 4B). KEGG pathway analysis revealed that these genes were significantly enriched in the PI3K-Akt signalling pathway, focal adhesion, and ECM-receptor interaction. Moreover, these genes were implicated in disease-associated pathways, including human papillomavirus infection, AGE-RAGE signalling in diabetic complications, and rheumatoid arthritis (Figure 4C). Collectively, these findings suggest that the identified genes may play critical roles in immune regulation, cell migration, and inflammatory responses.

PE-associated secreted proteins have been implicated in the pathogenesis and progression of DD. To identify potential biomarkers, an integrative analysis combining PE-associated secreted proteins with major depressive disorder (MDD)-related pathogenic genes was conducted, resulting in the identification of 15 candidate genes (Figure 5A). On the basis of these candidates, a diagnostic model was constructed to identify patients with comorbid PE and DD. LASSO regression was subsequently applied to refine the 15-gene set, yielding three potential pathogenic genes: CLEC3B, CTLA4, and PDPR (Figures 5B, C). In parallel, the RF machine learning algorithm was employed to rank candidate genes according to the mean decrease in the Gini index, and genes with MeanDecreaseGini ≥ 1.5 were selected (Figures 5D, E). An intersection analysis of the LASSO and RF results revealed three overlapping feature genes (Figure 5F).

Gene expression analysis of the GSE39653 dataset revealed significantly elevated levels of CLEC3B, CTLA4, and PDPR in DD patients than in healthy controls (Figure 5G) (P < 0.05). Similarly, in PE patients, the expression of these genes was also significantly upregulated relative to that in normal subjects (Figure 5H) (P < 0.05). These findings suggest that CLEC3B, CTLA4, and PDPR are consistently upregulated in both DD and PE patients and may serve as diagnostic biomarkers for identifying individuals with concurrent inflammatory or vascular dysfunction associated with depressive pathology.

To improve the diagnostic and predictive accuracy of DD associated with PE, a nomogram prediction model was developed on the basis of three selected feature genes via logistic regression analysis (Figure 6A). The cumulative scores corresponding to each feature gene were summed to calculate a total score, which was subsequently used to estimate the risk of PE-associated DD onset. Calibration curve analysis demonstrated high concordance between the predicted probabilities of the nomogram model and the ideal reference model (Figure 6B). Furthermore, decision curve analysis (DCA) was employed to evaluate the clinical utility of the model. The results indicated that the intervention strategy guided by the model yielded greater net benefit than did the intervention strategy in all patients or no intervention across most threshold probabilities (Figure 6C). The diagnostic sensitivity and specificity of the nomogram model for PE-associated DD were assessed via receiver operating characteristic (ROC) curve analysis, with an area under the curve (AUC) of 0.726 in the internal validation dataset GSE39653 (Figure 6D). In the external validation dataset GSE98793, the model achieved an AUC of 0.623 for diagnosing DD patients (Figure 6E), suggesting satisfactory predictive performance for depressive disorders related to inflammatory bowel disease.

GSEA was conducted to investigate the potential biological functions of three feature genes within the DD dataset GSE39653. The analysis revealed that the enriched KEGG pathways among the highly expressed feature genes predominantly involved the B-cell receptor, Toll-like receptor, and NOD-like receptor signalling pathways, along with antigen processing and presentation. The activation of the JAK-STAT and phosphatidylinositol signalling pathways further supported an immune-activated status, suggesting that these three genes may mediate immune regulation in DD (Figures 7A–C). The CIBERSORT algorithm was subsequently applied to estimate the infiltration levels of 22 immune cell types in DD samples. The results revealed a significant increase in naive CD4+ T cells (T_cells_CD4_naive) in the DD group (Figure 7D). Correlation analyses further demonstrated a negative association between CLEC3B and memory B cells (B_cells_memory), whereas PDPR was positively correlated with naive CD4+ T cells and plasma cells.

To elucidate the regulatory landscape of the three feature genes and their interacting components, a comprehensive gene–transcription factor (TF) regulatory network was constructed (Figure 8A). These feature genes are regulated by multiple transcription factors, suggesting that their functions in disease pathogenesis are tightly controlled by complex transcriptional programs. Accumulating evidence highlights the immunotoxic potential of plasticizers, underscoring their relevance to human immune modulation. To explore potential binding interactions between plasticizers and key genes associated with PE-related DD, a gene–compound interaction network was constructed via the Comparative Toxicogenomics Database (CTD) (Figure 8B). PDPR exhibited a stable binding affinity with bisphenol A (BPA), primarily forming hydrogen bonds at the ALA-3236 and GLY-239 residues (Figure 8C). CLEC3B interacted with both BPA and DEHP, with molecular docking analysis revealing minimum binding energies of –6.4 kcal/mol and –5.4 kcal/mol, respectively (Figures 8D, E). CTLA4 was found to interact with BPA and DBP. BPA formed hydrogen bonds at GLY-112 and ILE-109, whereas DBP bound at GLN-39 (Figure 8F, G). These findings indicate that these feature genes exhibit broad binding potential to various plasticizer components, indicating their toxicological involvement in the pathogenesis of PE-associated DD.