Bioinformatic Analysis of Apoptosis-Related Genes in Preeclampsia Using Public Transcriptomic and Single-Cell RNA Sequencing Datasets

Introduction

Preeclampsia (PE) is a pregnancy-specific complication mainly characterized by high blood pressure and proteinuria after 20 weeks of pregnancy. The global incidence of PE is approximately 2–8%, with higher rates in certain high-risk groups, such as women with first-time pregnancies, older age, multiple pregnancies, or a family history of PE.¹ This disease not only poses a serious threat to the health of pregnant women, but can also lead to intrauterine growth restriction of the fetus, iatrogenic premature delivery, and even endangers the life of the mother and fetus in severe cases.² The clinical manifestations of PE range from mild hypertension and proteinuria to severe renal dysfunction, liver dysfunction, retinopathy, and neurological complications, which vary in severity and speed of progression, and pose challenges for clinical diagnosis and management. Although previous studies have shown that PE may be related to genetics, immune system disorders, vascular endothelial dysfunction, and other factors, the specific pathogenesis of PE is not fully understood, there is still a lack of effective early prediction methods, and early termination of pregnancy is the only effective treatment is early termination of pregnancy.³ Therefore, early diagnosis and effective intervention based on pathogenesis are key to improving pregnancy outcomes in PE.

Apoptosis is a form of programmed cell death and is a widespread biological process. It encompasses related mechanisms such as phagocytosis, which is crucial for the clearance of apoptotic cells.⁴ Additionally, several other processes are closely related to apoptosis. Programmed necrosis (a regulated form of necrosis), pyroptosis (inflammation-related),⁵ autophagy (a self-degradation process that maintains intracellular homeostasis),⁶ and ferroptosis (an iron-dependent form of cell death)⁷ may all play significant roles in diseases such as preeclampsia.⁸ Apoptosis can help remove cells that are damaged or no longer required, thereby protecting the organism from potential pathological changes.⁹ Apoptosis plays a key role in the PE¹⁰ cells PE is a pregnancy complication that is closely related to abnormal placental function. The main placental cell types include villous cytotrophoblasts (VCT), syncytiotrophoblasts (SCT), and extravillous trophoblasts (EVT). These cells maintain normal placental functions through continuous differentiation and renewal. Studies have shown that the rate of placental cell apoptosis is abnormal in patients with PE, which may be related to insufficient blood flow and oxygen supply to the placenta.^11,12 In the context of PE, apoptosis of VCT and EVT cells may increase, while apoptosis of SCT cells may decrease, and this unbalanced apoptotic activity may lead to placental dysfunction, which in turn affects the entire course of pregnancy.^13,14

During a normal pregnancy, the apoptosis of placental cells is in a dynamic balance to maintain the normal development and function of the placenta. However, this balance is disrupted in patients with preeclampsia. Factors such as oxidative stress, inflammatory responses, and placental ischemia and hypoxia can induce abnormal activation of the apoptosis pathways in placental cells.¹⁵ For example, excessive reactive oxygen species (ROS) can damage intracellular DNA, proteins, and lipids, thereby activating endogenous apoptosis pathways, promoting the occurrence of caspase cascade reactions, and ultimately leading to cell apoptosis.¹⁶ Additionally, the abnormal elevation of inflammatory factors such as tumor necrosis factor-alpha (TNF-α) can activate death receptors on the cell surface through exogenous apoptosis pathways, triggering cell apoptosis.¹⁷ This abnormal apoptosis further disrupts the structure and function of the placenta, leading to abnormal placental vascular remodeling and decreased trophoblast invasion ability, ultimately resulting in a series of clinical symptoms of preeclampsia.¹⁷ Current research has identified various potential biomarkers associated with preeclampsia. In addition to the apoptosis-related genes (ARGs) focused on in this study, placental growth factor (PlGF) and soluble fms-like tyrosine kinase-1 (sFlt-1) have also garnered attention.¹⁸ PlGF is a cytokine that promotes angiogenesis, and its levels are typically significantly reduced in patients with preeclampsia, which is closely related to placental vascular dysplasia.¹⁹ sFlt-1 is an anti-angiogenic factor that is highly expressed in patients with preeclampsia; it can bind to angiogenic factors like PlGF, blocking their biological activity and leading to an imbalance in angiogenesis, which is one of the important molecular mechanisms of preeclampsia.²⁰ In terms of molecular regulatory mechanisms, there exists a complex interaction network among these biomarkers. For instance, certain genes within ARGs may influence the expression of PlGF and sFlt-1 by regulating intracellular signaling pathways, thereby affecting placental vascular generation and function. Currently, commonly used clinical diagnostic methods include measuring the levels of related biomarkers in the blood of pregnant women and ultrasound examinations to assess the condition of the placenta and fetus.²¹ However, these methods have certain limitations; therefore, in-depth exploration of new diagnostic biomarkers and molecular regulatory mechanisms is of great significance for improving the accuracy of early diagnosis of preeclampsia.

Although apoptosis plays an important role in the occurrence and development of PE, the specific molecular mechanism remains unclear, prompting researchers to explore the specific biological role of apoptosis in PE and its molecular mechanisms with a view to finding possible therapeutic targets to improve pregnancy outcomes. In summary, the study of apoptosis in PE not only contributes to our understanding of the nature of this complex pathological state but may also provide key information for the development of new prevention and treatment strategies.

Recent studies have identified several biomarkers associated with preeclampsia. For example, CD93, a marker related to phagocytosis, has been shown to be associated with preeclampsia.²² CD93 is involved in various cell signaling pathways. It can interact with various ligands and receptors on the cell surface, affecting processes such as cell adhesion and migration.²³ In the context of apoptosis and immune response, CD93 may regulate the clearance of apoptotic cells by immune cells, and its abnormal expression in preeclampsia may disrupt the normal balance of placental cell renewal and immune regulation, thereby promoting the development of the disease.²⁴ In addition to CD93, other potential biomarkers for the diagnosis of preeclampsia have been studied, such as TRA4, PROCR, and TXNIP, which were identified by Li Weiwen et al through bioinformatics and machine learning methods as potential PE biomarkers.²⁵ These continuously discovered biomarkers deepen our understanding of the pathogenesis of preeclampsia. Furthermore, apoptosis is closely related to oxidative stress, which is also associated with the development of preeclampsia.²⁶ Xanthine oxidase is an enzyme involved in the generation of reactive oxygen species (ROS) during oxidative stress.²⁷ The serum uric acid to serum creatinine ratio, as an indirect marker of xanthine oxidase activity, has been shown to be related to the development of preeclampsia.²⁸ Elevated levels of xanthine oxidase and abnormal serum uric acid to serum creatinine ratios may lead to increased oxidative stress in the placenta, triggering placental cell apoptosis, disrupting normal placental function, and ultimately promoting the occurrence and development of preeclampsia.²⁹ Thus, research on apoptosis in preeclampsia not only helps us understand the nature of this complex pathological state but may also provide key information for developing new prevention and treatment strategies.

Based on transcriptome and single-cell RNA sequencing (scRNA-seq) data in public databases, this study explored the potential diagnostic biomarkers and molecular regulatory mechanisms of PE through apoptosis-related genes using bioinformatics to provide new ideas for the clinical diagnosis, prevention, and treatment of PE.

Materials and Methods

Data Source

Two transcriptome datasets related to PE and a single-cell RNA sequencing (scRNA-seq) dataset were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds). The GSE43942 dataset, based on the GPL10191 platform using a chip sequencing method, served as the training set for this study and included placental tissue samples from 5 PE patients and 7 normal individuals.³⁰ The clinical information table of the GSE43942 dataset was shown in the Table S1. The GSE25906 dataset, utilized as the validation set, was based on the GPL6102 platform with chip sequencing, featuring placental tissue samples from 23 patients with PE and 37 normal individuals.³¹ The clinical information table of the GSE25906 dataset was shown in the Table S2. The scRNA-seq dataset GSE173193, based on the GPL24676 platform with high-throughput sequencing, included placental tissue samples from 2 PE patients and 2 normal individuals.^32,33 The clinical information table of the GSE173193 dataset was shown in the Table S3. In addition, 680 apoptosis-related genes (ARGs) were extracted for analysis from a previously published article for analysis.³⁴

ScRNA-Seq Analysis

To enhance our understanding of PE at the cellular level, we used the Seurat package (v 5.0.1)³⁵ to perform a comprehensive analysis of the scRNA-seq data in the GSE173193 dataset. Strict quality control (QC) measures were implemented on the scRNA-seq data to exclude low-quality data that might have been caused by cell damage or library preparation defects. To filter the single-cell sequencing data, the data was first integrated. GSM5261695 and GSM5261696 were used as the control group, and GSM5261699 and GSM5261700 were used as the disease group. After being read by read10X, the data was converted into a Seurat object. During the reading process, cells with fewer than 100 genes and genes covered by fewer than 3 cells were filtered out. The PercentageFeatureSet function in the “Seurat” R package was used to calculate the “nCount/nFeature-RNA” and “percent.Mt” of the samples for data quality control. The “Mt” of the samples was calculated by the percentage feature set function. Next, the percentages of the number of genes (nfeature), the number of cells (ncount), and the number of mitochondrial sequencing reads (Mt) in each sample were calculated. Meanwhile, cells with more than 20% mitochondrial genes were filtered out. Cells with the number of intracellular genes less than or equal to 100 and greater than or equal to 10,000 were removed, and genes with count numbers less than or equal to 100 and greater than or equal to 30,000 were also removed, so as to exclude the cases of multi-cell bodies or some cells with poor quality or in a dying state. The 24,100 cells and 22,891 genes that met the screening criteria were used for subsequent analysis. Subsequently, the data were normalized using the NormalizeData function, and the top 2500 highly variable genes were selected for downstream analysis using the FindVariableFeatures function based on the relationship between mean and variance. The results were visualized by applying the LabelPoints function, which highlighted the top ten genes with the greatest variation. The data were further normalized using the ScaleData function, and an elbow plot was created to determine the appropriate principal components (PCs) for further analysis. Utilizing the standard Seurat workflow, Uniform Manifold Approximation and Projection (UMAP) was employed to perform clustering analysis on dimensionality-reduced cells, with the resolution set at 0.5. After obtaining cell clusters through clustering, marker genes were identified from the GSE173193 dataset according to the reference¹² used in this study. Based on these markers, cell annotation was conducted by classifying cells that exhibited high specificity for gene expression within clusters. For clusters lacking marker genes, annotations were supplemented using singleR (v 2.0.0),³⁶ with the threshold set at a quantile of 0.7. Under these circumstances, cell types specific to this study were annotated and a histogram was created to display the proportions of different cell types in the PE and control samples.

Differential Expression Analysis

When conducting differential expression analysis, since a large number of genes were simultaneously tested, in order to avoid an excessive number of false-positive results due to multiple testing, we adopted a strict statistical correction method. Specifically, we used the Benjamini-Hochberg method to control the False Discovery Rate (FDR).³⁷ This method sorted all the P-values of the tests and calculated the corrected q-values based on the sorting results to ensure that the proportion of false discoveries during the entire testing process was maintained at a reasonable level. In actual operation, we took the corrected q-value as one of the important bases for judging the differential expression of genes. Only when the q-value was less than the pre-set threshold (such as 0.05) was the differential expression of the gene considered to be statistically significant.³⁸ Subsequently, differential expression analysis was performed. Initially, within the GSE173193 dataset, differences between the PE and control samples across various cell types were compared using the Wilcoxon test, and box plots were used for visualization. Cell types showing significant differences in abundance between samples were defined as differential cells for subsequent analysis (P < 0.05). In addition, in the GSE173193 dataset, the FindMarkers function was used to screen for differentially expressed genes (DEGs) within differential cells between PE and control samples, with thresholds set at |log₂Fold Change (FC)| ≥ 0.5, adj. P < 0.05, and ptc > 0.1, defined as DEGs1, and a Manhattan plot was created for visualization. In the training set, differential expression analysis was performed using the limma package (v 3.54.0)³⁹ to identify DEGs between PE and control samples, with the thresholds |log₂FC| > 0.5 and P < 0.05, defined as DEGs2. For more intuitive observation of gene differences, volcano plot and heatmap were generated using ggplot2 (v 3.4.1)⁴⁰ and ComplexHeatmap (v 2.15.1),⁴¹ respectively, where the volcano plot displayed the top 10 up- and down-regulated genes sorted by log₂FC.

Identification and Enrichment Analysis of Differentially Expressed ARGs (DE-ARGs)

Previously identified DEGs1, DEGs2, and ARGs were intersected using the VennDiagram package (v 1.7.1),⁴² resulting in the identification of DE ARGs. These DE-ARGs were then subjected to Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses to explore their biological functions and signaling pathways. This process was performed using the clusterProfiler package (v. 4.2.2),⁴³ with a significance level of P < 0.05. The GO terms included three categories: biological processes (BP), molecular functions (MF), and cellular components (CC). To explore the protein-level interactions of DE-ARGs, these genes were input into the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org) with a confidence score threshold of > 0.15 for selection. The resulting protein-protein interaction (PPI) network was visualized using the Cytoscape software (v. 3.5.2).⁴⁴

Recognition and Expression Assessment of Biomarkers

In the training set, using the constructed PPI network, five algorithms from the CytoHubba plugin in Cytoscape software were employed: degree, MCC (Maximal Clique Centrality), DMNC (Density of Maximum Neighborhood Component), MNC (Maximum Neighborhood Component), and EPC (edge-percolated component) were employed. Each algorithm was used to identify the top 9 genes based on their respective scores. The intersection of these genes across all the algorithms yielded a set of hub genes, and Venn diagrams were generated using the VennDiagram package for visual representation. The expression of hub genes in the PE and control samples was verified using both training and validation datasets. Genes exhibiting consistent expression trends and significant differences between the two datasets were defined as biomarkers for subsequent analyses (P < 0.05).

Construction of a Nomogram Through Biomarkers

To maximize the diagnostic value of biomarkers for PE, the rms package (v. 6.5–1)⁴⁵ was used to construct the nomogram. Within this nomogram, a scoring system was devised based on the expression levels of biomarkers in the training set, with each factor assigned a specific score. The cumulative sum of the scores for various factors represents the total number of points. Based on the total points, the probability of a PE diagnosis was predicted, with higher scores indicating a greater likelihood of PE. Furthermore, a calibration curve was plotted to perform the Hosmer-Lemeshow (HL) test, a critical diagnostic for the predictive accuracy of the nomogram, and a p-value greater than 0.05, from the HL test, indicated that the predicted values did not significantly deviate from the actual values, thus passing the test. The ggDCA package (v 1.2) (https://www.rdocumentation.org/packages/ggDCA/versions/1.1) was used to conduct decision curve analysis (DCA) to assess the clinical utility of the nomogram.

Gene Set Enrichment Analysis (GSEA) of Biomarkers

After identifying the biomarkers, we conducted a series of analyses to characterize their functionality. To uncover the biological functions and pathways linked to biomarkers, the “h.hallmark.v7.4. symbols.gmt” from the Gene Set Enrichment Analysis (GSEA) database (http://www.broadinstitute.org/gsea/index.jsp) was selected as the reference gene set. Spearman correlation coefficients were calculated between each biomarker and all other genes in the training set, ranking the genes in descending order based on these coefficients to identify the related genes for each biomarker. Subsequent GSEA for each biomarker, conducted via the clusterProfiler package, considered results with P < 0.05 as significantly enriched. The top five enriched pathways were visualized using the enrichment plot package (v. 1.18.3).⁴⁶

Functional and Annotation Analysis

Variations caused by genetic and non-genetic diseases can be explored through the association between genotypes and phenotypes.⁴⁷ The human genome contains hundreds of single-nucleotide variations (SNVs), many of which are known to cause diseases.⁴⁸ Approximately 21% of the amino acid variations lead to missense mutations at different protein sites, known as post-translational modifications, which can result in diseases.⁴⁹ Consequently, to understand the SNV profiles of biomarkers, the ActiveDriverDB online tool (https://activedriverdb.org/) was used for SNV analysis of biomarkers.

Additionally, analysis was performed to evaluate the functional similarity of biomarkers using the GOSemSim package (v.2.27.2).⁵⁰ The required GO annotation data for human species across the BP, CC, and MF aspects were constructed using GOSemSim. Semantic similarities were calculated using the mgeneSim function. Thereafter, the geometric mean of these similarities across all categories was computed to derive a final score to evaluate the overall functional impact of the biomarkers across the three GO domains.

GeneMANIA (http://www.genemania.org/) was used to predict genes functionally related to biomarkers, and a co-expression network of the top 20 interacting genes and their associated functions was visualized.

Physicochemical properties of proteins reflect the fundamental characteristics of an organism.⁵¹ To understand the physicochemical properties of the biomarkers, the amino acid sequences were retrieved from UniProt (https://www.uniprot.org/). Basic information on the protein-coding genes (biomarkers), including amino acid composition, molecular weight, protein half-life, and calculated parameters, such as the protein instability index and theoretical isoelectric point, was analyzed using the ProtParam tool (https://web.expasy.org/protparam/).

Immune Infiltration and Pathway Analyses

Immune infiltration assays were performed to explore variations in the immune environment between patients with PE and healthy controls. Specifically, the CIBERSORT algorithm (v 1.03)⁵² was used to estimate the abundance of 22 immune cell types³⁰ in each sample of the training set. Samples were filtered based on reliability, excluding those with zero expression and removing samples with P > 0.05. Further analysis was performed to examine the distribution of these immune cells in PE versus control samples, with a significance threshold of P < 0.05. Spearman correlation analysis was conducted using the psych package (v 2.2.9)⁵³ to identify potential relationships between significantly different immune cells and biomarkers (|cor| > 0.3 and P < 0.05).

In addition, 13 immune-related pathways previously identified in the literature⁵⁴ were explored to assess the differences in pathway activity between patients with PE and control samples in the training set. The single-sample GSEA (ssGSEA) algorithm within the GSVA package (v 1.42.0)⁵⁵ was used to calculate single-sample scores for these immune-related pathways, allowing for a comparison of immune-related pathway scores between the PE and control samples (P < 0.05). Further analysis was conducted using the psych package to perform Spearman correlation analysis between these differentially immune-related pathways and biomarkers (|cor| > 0.3 and P < 0.05).

Construction of Regulatory Network

Further analyses were performed to elucidate the molecular regulatory mechanisms that affect these biomarkers. Using the multiMiR package (v 1.18.0),⁵⁶ miRNAs associated with biomarkers were predicted using default settings (predicted.cutoff and predicted.cutoff.type=“p”, selecting the top 20%). Key miRNAs were consistently predicted across five databases: DIANA-microT (http://www.microrna.gr/microT), ElMMo (http://www.mirz.unibas.ch/ElMMo3/), miRecords (http://c1.accurascience.com/miRecords/), TargetScan (http://www.targetscan.org/) and miRDB (http://www.mirdb.org/). Subsequent predictions of potential lncRNA interactions with these miRNAs were performed using the starBase (https://starbase.sysu.edu.cn/) and miRNet (https://www.mirnet.ca/) platforms. In starBase, a threshold of pancancerNum > 4 was used to filter the data, retaining only those lncRNAs that were consistently predicted by both databases to have interactive relationships. Using this information, an lncRNA-miRNA-mRNA regulatory network was constructed using the Cytoscape software. Similarly, the ENCODE database (https://www.encodeproject.org/) within miRNet was used to identify transcription factor (TFs)-targeting biomarkers and a TF-mRNA network was constructed using Cytoscape software.

Drug Prediction and Molecular Docking Analyses

In this study, potential drug-targeting biomarkers were identified using the enrichR package (v. 3.1)⁵⁷ from the Drug SIGnatures database (DSigDB) (http://dsigdb.tanlab.org/DSigDBv1.0/). Results with an adj.P < 0.05, were considered statistically significant and potential drugs, and the top 10 drugs, based on scoring, were selected for visualization. In addition, the relationship between differential immune cells obtained from immune infiltration analysis and biomarkers was incorporated into an interaction network of drugs, biomarkers, and immune cells using Cytoscape software, which facilitated the understanding of the interactions among drugs, immune cells, and biomarkers.

To understand the stereochemical properties of the biomarkers, AlphaFold v2.0 (https://alphafold.ebi.ac.uk/) was employed as an artificial intelligence tool to predict the 3D structures of the biomarkers. The amino acid sequences of these biomarkers were downloaded from the UniProt-KB database (http://www.uniprot.org/) and the local distance difference test (LDDT) scores from the AlphaFold database were used to assess the stereochemical performance of these 3D models. To further explore the role of biomarkers in drug therapy, the molecular structures of drugs and protein 3D structures were downloaded from PubChem (https://pubchem.ncbi.nlm.nih.gov/) and Protein Data Bank (PDB) (http://www.rcsb.org/) databases, respectively. CB-Dock2 was used for the molecular docking of drugs against biomarkers to obtain the binding energies. It is generally accepted that the more stable the ligand-receptor conformation, the lower the required binding energy, with energies ≤ −5 kcal/mol considered indicative of strong binding affinity. Drugs with the highest 2D conformation scores were selected for molecular docking, with visualization conducted using PyMOL software (v. 3.1.1).⁵⁸

Cell Enrichment and Communication Analyses

Among the previously annotated cell types, a series of analyses were conducted to explore the cellular mechanisms of PE across all samples in GSE173193. Applying the ReactomeGSA package (v 1.12.0),⁵⁹ functional enrichment was performed on each of the identified differential cells from the prior differential analysis to explore their biological pathways, with the top 20 pathways displayed. Subsequently, with support from the celltalker package (v 0.0.7.9000),⁶⁰ PE and control samples were distinguished within the training set to conduct cell communication analysis. This involves detecting the expression and pairing of receptors and ligands within these cell types, with the aim of inferring interactions between different cells.

Identification of Key Cell and Pseudo-Time Analysis

Biomarkers were examined for their expression differences in differential cells between PE and control samples. Violin plots were used for visualization, identifying cells with significant expression differences in most biomarkers as key cells. Further analysis illustrated the distribution of biomarkers across the different cell types.

To investigate the developmental trajectory and evolution of key cells, the cells were initially subjected to dimensionality reduction and clustering to identify cellular subtypes. Subsequently, the monocle package (v 2.26.0)⁶¹ was employed to conduct a pseudo-time series analysis using 2500 highly variable genes across both the PE and control samples. In addition, temporal changes in biomarker expression within key cells were monitored.

Reverse Transcription-Quantitative Polymerase Chain Reaction (RT-qPCR)

Placental tissues from six patients (three in the control group and three in the preeclampsia group) were collected after delivery at the First People’s Hospital of Yunnan Province. Informed consent was obtained from all participants. This study was approved by the Ethics Committee of First People’s Hospital of Yunnan Province (approval number: 221). The RNA extracted from each sample was quantified using 1 mL of triazole reagent and 300 µL of chloroform, the concentration of which was quantified using a nanophotometer N50. cDNA was reverse-transcribed using the first-strand cDNA of SweScript. Synthetic Kit (Service Biology, Wuhan, China). RT-qPCR amplification was performed for 40 cycles using a CFX96 real-time fluorescent quantitative PCR instrument (Bio-Rad, Hercules, California, USA). The PCR cycle consisted of 1 min of pre-denaturation at 95°C, 20s of denaturation at 95°C, 20s of annealing at 58°C, and 30s of extension at 72°C. During RT-qPCR analysis, each reaction mixture contained 3 µL CDN A,5 microliters µL 2x Universal Blue SYBR GreengPCR master mixture, 1 µL forward primer, and 1 µL reverse primer. The sequence information for the four gene primers is in Table 1.

Table 1 Sequence Information Table for the Four Gene Primers

Statistical Analysis

All analyses were performed using R software (v 4.2). Differences between the groups were analyzed using the Wilcoxon test. Statistical significance was defined as P < 0.05, with *P < 0.05, **P < 0.01, ***P < 0.001, ****P < 0.0001, indicating increasing levels of significance, and ns denotes non-significant (P > 0.05).

Results

Annotation of 11 highly Specific Cell Types

For the GSE173193 dataset, the results before and after quality control are illustrated in Figure S1A and B, revealing that 24,100 cells and 22,891 genes met the criteria. From these, 2500 highly variable genes were selected for downstream analysis (Figure S1C). Subsequently, stability was achieved after the first 20 PCs; hence, the top 20 PCs were selected for further analysis (Figure S1D and E). Based on this, UMAP clustering identified 21 cell clusters, and their distribution in PE and control samples is displayed (Figure 1A and B). The bubble plot showed a high specificity for marker genes, leading to the cells being named based on these markers (Figure 1C). Ultimately, 11 cell types were annotated: granulocytes, villous cytotrophoblasts (VCT), macrophages, myelocytes, extravillous trophoblasts (EVT), T/NK cells, monocytes, syncytiotrophoblasts (SCT), common myeloid progenitors (CMP), endothelial cells, and erythroblasts (Figure 1D). Histograms displayed the proportions of these cell types in PE and control samples, clearly showing that VCT was more prevalent in PE samples, whereas granulocytes were more abundant in the control samples (Figure 1E).

Figure 1 Annotation maps for 11 highly specific cell types. (A) 21 samples UMAP cell taxon clustering, Different colors represent different cell populations. (B) PE and control grouping UMAP cell taxon clustering. (C) Bubble plots of marker genes in different cellular taxa, darker colors represent higher gene expression. (D) Clustering map of 11 maker genes. (E) Histogram of the percentage of 11 cell types between PE and control groups: different colors represent different cell types.

Recognition and Functional Characterization of DE-ARGs

In the GSE173193 dataset, differences between the PE and control samples were compared across 11 cell types, revealing significant variations in 10 cell types (P < 0.05), excluding macrophages (P > 0.05), which were defined as differential cells (Figure 2A). Further analysis of 10 differential cells identified 5577 DEGs between PE and control samples. After removing duplicates, 3674 unique DEGs1 were identified, with 790 genes upregulated and 2884 downregulated in the PE samples (Figure 2B). Similarly, in the transcriptome dataset (training set), 1882 DEGs2 were identified between PE and control samples, with 810 genes upregulated and 1072 downregulated in the PE samples (Figure 2C and D). Subsequently, an intersection of 3674 DEGs1, 1882 DEGs2, and 680 ARGs was obtained, yielding 28 DE-ARGs (Figure 2E), which were further included in the enrichment analysis. This analysis identified 779 GO entries, comprising 673 BP, 26 CC, and 80 MF, including “regulation of the apoptotic signaling pathway”, “intrinsic apoptotic signaling pathway”, “extrinsic apoptotic signaling pathway via death domain receptors”, “negative regulation of the apoptotic signaling pathway”, “intrinsic apoptotic signaling pathway in response to DNA damage”, “NAD+ nucleosidase activity”, and “vesicle lumen” (Figure 2F). Additionally, 56 KEGG pathways were identified, such as the “AGE-RAGE signaling pathway in diabetic complications”, “apoptosis”, “chagas disease”, “focal adhesion”, and “lipid and atherosclerosis” (Figure 2G). These enrichment results, involving cell death, signaling pathway regulation, and extracellular interactions, revealed the potential processes through which DE ARGs affected PE.

Figure 2 Identification and Functional Characterization Map of DE-ARGs. (A) 11 Differential boxplots of cell types. “ns” represented not significant, “*” represented p < 0.05, “***” represented p < 0.001, and “****” represented p < 0.0001. (B) Differential Gene Manhattan Map: Red dots represent up-regulated genes and green dots represent down-regulated genes; the 5 genes with the most significant difference between up- and down-regulation show the gene names. (C) Volcano plot of gene expression in PE and control samples, red represents high expression, blue represents low expression, and the shade of the color reflects the amount of expression. (D) Heat map distribution of PE and control samples, the top section is a heatmap of the expression density of up- and down-regulated differential genes in the samples, showing the lines of the five quartiles and the mean; the bottom section is a heatmap of the expression of the differential genes. (E) Wayne diagrams between DEG1, DEG2 and ARGs. (F) GO enrichment results: the horizontal axis is the number of genes enriched to any one pathway, and the vertical axis is the name of the GO-enriched entry. (G) KEGG enrichment results: gray dots indicate genes, yellow dots indicate different pathway names, and different colored lines represent pathways that genes are enriched into.

Incorporating Biomarkers CD44, MIF, PIK3R1, and TLR4 into a Nomogram

Based on the 28 DE-ARGs, a PPI network with 22 nodes and 38 edges was constructed, where JUN, TLR4, CD44, and ICAM1 had the highest degrees (Figure 3A). Further analysis was conducted using 5 different algorithms from the CytoHubba plugin. As DMNC only scored the top nine genes, we selected the top nine genes from each algorithm and identified their intersection, resulting in nine hub genes: CD44, ICAM1, JUN, MIF, PARP1, PIK3R1, SERPINE1, THBS1, and TLR4 (Figure 3B). These nine hub genes were included in the expression level analyses, which clearly revealed consistent expression trends and significant intergroup differences in both the training and validation datasets for CD44, MIF, PIK3R1, and TLR4 (P < 0.05) (Figure 3C and D). Notably, CD44 and TLR4 were expressed at significantly lower levels in the PE samples (P < 0.05), whereas MIF and PIK3R1 were expressed at higher levels (P < 0.05). These four genes were identified as potential biomarkers in this study.

Figure 3 Selection, expression and validation of biomarkers. (A) PPI network of 28 DE-ARGs, the size of the dots represents the number of cells of that type, and the thickness of the lines represents the strength of communication between the corresponding cell groups. (B) 5 Algorithms Top10 Gene Intersection. (C) Boxplot of biomarkers expression in GSE43942. The blue and red represented the “control” group and the “PE” group, respectively. “*” represented p < 0.05, and “**” represented p < 0.01. (D) Boxplot of biomarkers expression in GSE25906. “ns” represented not significant, “**” represented p < 0.01, and “***” represented p < 0.001. (E) The column-line diagram model of Biomarker. (F) Calibration curves for line diagram models: the horizontal coordinates of the calibration curve for the column-line diagram model indicate the probability of illness predicted by the nomogram, and the vertical coordinates indicate the actual probability of illness. (G) DCA Decision Curve: the horizontal coordinate is the threshold probability and the vertical coordinate is the net benefit rate after benefits minus drawbacks.

Subsequent analysis integrated these biomarkers into a nomogram to maximize their diagnostic value for PE (Figure 3E). The calibration curve compared the predicted probabilities of the nomogram with the observed actual probabilities (Figure 3F). The P-value of the HL test was 0.489, indicating no significant difference between the predicted and actual values, and the mean absolute error (MAE) was 0.037, suggesting a minimal error between the actual and predicted risk of disease. This demonstrates the high accuracy of the nomogram in predicting the risk of illness in the sample. DCA also indicated that the nomogram’s net benefit was higher than those of the positive and negative controls, providing greater net benefits than the use of any single biomarker alone, thereby exhibiting substantial clinical utility (Figure 3G).

Investigating the Potential Functions and Involved Signaling Pathways of Biomarkers

GSEA was performed on biomarkers to explore enriched signaling pathways. The results revealed that CD44, MIF, PIK3R1, and TLR4 were enriched in 24, 14, 30, and 10 hallmark pathways, respectively (Figure 4A–D). In particular, they were commonly enriched in the “epithelial mesenchymal transition” pathway. CD44, MIF, and PIK3R1 were jointly enriched in “E2F targets”; CD44, PIK3R1, and TLR4 were commonly found in “MYC targets V1”; and CD44 and TLR4 were both enriched in “oxidative phosphorylation”. These pathways are closely related to the regulation of cellular functions and pathological states, suggesting the mechanisms by which these biomarkers might influence the development and progression of PE.

Figure 4 Potential functions of biomarkers and associated signaling pathways. (A) CD44 enrichment analysis, (B) MIF enrichment analysis, (C) PIK3R1 enrichment analysis, (D) TLR4 enrichment analysis. Each fold represents a pathway with lines marking the genes located in the gene set. Analysis of (E) CD-44, (F) MIF, (G) PIK3R1, (H) TLR4 single nucleotide variants. (I) Functional similarity analysis of biomarkers. (J) Biomarker GGI Network: different colored lines represent different interactions, and different color blocks represent the different functional roles involved.

Further research on SNVs of biomarkers revealed significant findings. In CD44, 65.9% of the sequences were predicted to be disordered, indicating that these regions lacked a fixed three-dimensional structure and might have played a crucial role in regulating protein function and protein-protein interactions (Figure 4E). The mutations affected various aspects, including direct effects, network rewiring, motif changes, and proximal and distal associations, which could have altered how the protein interacted with other molecules, thereby changing signal transduction pathways and cellular functions. The sites were related to acetylation, methylation, succinylation, SUMOylation, ubiquitination, glycosylation, and phosphorylation, all of which could affect the activity, stability, interaction partners, and cellular localization of the protein. Similar results were observed for other biomarkers, with 0.0%, 10.08%, and 5.48% of the sequences in MIF, PIK3R1, and TLR4, respectively, predicted to be disordered, reflecting potential differences in their cellular functions and stabilities (Figure 4F–H). These results suggest that the sequences and functions of these biomarkers are regulated by various molecular mechanisms and that SNVs can significantly alter their functions.

Friends analysis evaluated the functional similarity among biomarkers, indicating that the four genes shared certain similarities, with TLR4 showing the highest overall functional similarity (scores > 0.5) (Figure 4I). GeneMANIA provides a co-expression network for these biomarkers, in which physical interactions are the most prevalent (77.64%), followed by co-expression effects (8.01%). These genes were collectively involved in functions such as “response to molecule of bacterial origin”, “cellular response to molecule of bacterial origin”, “positive regulation of tumor necrosis factor superfamily cytokine production”, “positive regulation of cell activation”, “regulation of phosphatidylinositol 3-kinase signaling”, “cellular response to biotic stimulus”, and “positive regulation of leukocyte activation” (Figure 4J). The physicochemical properties of these biomarkers are listed (Table 2). For example, CD44 has the formula C1719H2701N487O563S18 with a molecular weight of 39,775.36, consisting of 362 amino acids, a theoretical isoelectric point (pI) of 5.38, an aliphatic index of 72.46, a grand average of hydropathicity (GRAVY) of 0.449, an estimated protein half-life of 30 h, and an instability index of 32.11.

Table 2 Physical and Chemical Properties of the Biomarkers

Differential Immune Cell Infiltration and Pathway Activation in PE

Subsequent analysis aimed to explore the differences in the immune microenvironment composition between patients with PE and healthy individuals, initially showing the infiltration abundance of 22 immune cells in both types of samples within the training set (Figure 5A), with box plots further visualizing these disparities (Figure 5B). The results revealed significant differences in the two immune cells between the PE and control samples (P < 0.05). Monocytes had significantly less infiltration in PE patients (P < 0.05), while resting NK cells exhibited a higher infiltration abundance (P < 0.05). In addition, biomarkers showed varying degrees of correlation with these differentially infiltrated immune cells, with PIK3R1 displaying the highest negative correlation with monocytes (cor = −0.615 and P< 0.05) and TLR4 showing the highest positive correlation with resting NK cells (cor = 0.650, P< 0.05) (Figure 5C). These results suggest differences in the composition of immune cells within the immune microenvironment between PE patients and healthy individuals. The two identified immune cells might play roles in the pathogenesis of PE, and the identified biomarkers could potentially represent the immune microenvironment of patients with PE, which is crucial for their management.

Figure 5 Differential immune cell infiltration and pathway activation in PE. (A) Map of immune-infiltrating cell abundance, different colors represent different immune cells. (B) Box plot of immune cell differences between PE and control groups, horizontal coordinates indicate immune cells and vertical coordinates indicate immune cell immunity scores in the samples; red indicates the control group and blue indicates the PE control group. “ns” represented not significant, and “*” represented p < 0.05. (C) Heatmap of Differential Immune Cell-Biomarker Correlation: * denotes significance, numbers represent correlation coefficients, absolute value of correlation coefficients <0.3 is null. “*” represented p < 0.05, and “****” represented p < 0.0001. (D) Comparative box plots of immune pathway differences. “ns” represented not significant, “*” represented p < 0.05, and “**” represented p < 0.01. (E) Heatmap of biomarkers correlating with immune pathways. “ns” represented not significant, “*” represented p < 0.05, “**” represented p < 0.01, and “***” represented p < 0.001.

After completing the immune infiltration analysis, we further explored the differences in 13 immune-related pathways between patients with PE and healthy individuals in the training set. This has allowed for a deeper understanding of the pathogenesis of PE, revealing how these cells interact at the molecular level, affect disease progression, and identify potential therapeutic targets. The results showed that all four significantly different immune-related pathways, APC co-stimulation, checkpoint, T cell co-stimulation, and MHC class I, had notably lower scores in the PE samples (P < 0.05) (Figure 5D). This suggests that immune regulation in patients might be suppressed in some key pathways that are typically associated with immune activation and antigen presentation. The biomarkers also showed strong correlations with these pathways; for example, the correlation between MIF and checkpoints was as high as −0.804 (P < 0.01), and TLR4 had the strongest positive correlation with MHC class I (cor = 0.867 and P < 0.001) (Figure 5E). The reduced activity of these pathways could be a key factor in the development of PE, and the association between these pathways and biomarkers could provide potential targets for the development of targeted therapies for PE.

Revealing the Potential Molecular Regulatory Mechanisms of Biomarkers

Four miRNAs consistently predicted across five databases (DIANA-microT, ElMMo, miRecords, TargetScan, and miRDB) were identified as key miRNAs, whereas 16 lncRNAs consistently predicted in two other databases (starBase and miRNet) were identified as key lncRNAs. Based on this, an lncRNA-miRNA-mRNA regulatory network was constructed consisting of 21 nodes and 45 edges (Figure 6A). Therefore, PIK3R1 was the only mRNA with intersecting miRNA results, and lncRNAs such as NEAT1, OIP5-AS1, GAS5, and DNAAF4-CCPG1 may simultaneously regulate PIK3R1 through hsa-mir-29a-3p, hsa-mir-29b-3p, hsa-mir-29c-3p, and hsa-mir-221-3p.

Figure 6 The potential molecular regulatory mechanisms of biomarkers. (A) Biomarker-mRNA-lncRNA network: red is biomarker, green is miRNA and blue is lncRNA. (B) Biomarker-TF Network, red is biomarkers, blue is TF. (C) TOP10 Drug prediction network: biomarkers in red, Top 10 drugs in blue. (D) Drug-Biomarker-Immune cell prediction networks: green is an immune cell with differences. (E) TLR4 and (F) CD44 Molecular Docking: the left image is global, the right image is localized.

We constructed a TF mRNA network to further elucidate the potential regulatory mechanisms of these biomarkers. In this network, 67 TFs were predicted to be associated with three biomarkers, comprising 70 nodes and 75 edges (Figure 6B). TFs such as MXD3 and ZNF71 may co-regulate MIF and PIK3R1, whereas PHF8, SAP30, KDM5B, CBFB, and WRNIP1 may co-regulate MIF and CD44, and TARDBP may co-regulate CD44 and PIK3R1.

Exploring the Structure and Function of Drug Targets Associated with Biomarkers

In DsigDB, 312 potential drug targets identified from biomarkers were visualized, selecting the top 10 ranked drugs, including arsenic, fluorescein isothiocyanate (FITC), isoflavone, heparitin, Healing, acetaldehyde, indomethacin sodium, nandrolone phenpropionate, lorazepam, and 9001-31-4. A network comprising 14 nodes and 22 edges was constructed with ARSENIC targeting 4r biomarkers (Figure 6C). Moreover, the previously identified differential immune cells (monocytes and resting NK cells) were integrated to construct a network comprising 16 nodes and 26 edges, featuring interaction pairs such as isoflavone-TLR4-monocytes (Figure 6D). This suggests that isoflavone may influence monocyte function through TLR4, potentially affecting immune regulation and inflammatory responses. Subsequent analysis revealed structural models of the proteins corresponding to the biomarkers (Figure S2A–D). Notably, the overall color of the MIF model was blue, with a predicted LDDT score greater than 90, indicating that its predicted structure closely resembled the actual structure. Furthermore, the highest-scoring 2D conformations among the aforementioned drugs were selected for molecular docking, and both CD44 and TLR4 were docked using FITC. Specifically, the binding energy between CD44 and FITC was −8.1 kcal/mol, with residues such as ARG-29, SER-58, PHE-56, THR-133, and ASP-134 forming hydrogen bonds with FITC (Figure 6E). Similarly, the binding energy between TLR4 and FITC was −8.0 kcal/mol, with residues SER-126 and GLY-124 forming hydrogen bonds with FITC (Figure 6F).

Probing Cellular Communication and Biological Functions Under PE Conditions

Following transcriptomic analysis, an understanding of the expression patterns and functions of the biomarkers was obtained. To further explore the heterogeneity of these processes at the cellular level, biomarker expression profiles were analyzed from a single-cell perspective to explore unique expression patterns across different cell types. This approach aimed to further understand the mechanisms of PE at the cellular level using the scRNA-seq dataset (GSE173193).

Initially, enrichment analysis of ten different cell types revealed the main biological functions of these cells. These cells were enriched in pathways such as “regulation of thyroid hormone activity”, “hydroxycarboxylic acid-binding receptors”, “mitochondrial uncoupling”, “proton-coupled neutral amino acids”, and “classical antibody-mediated complement activation” (Figure S3A–J). These pathways play crucial roles in maintaining internal homeostasis and in responding to various physiological and pathological conditions.

Further analysis of cellular communication, displaying the number and weight of interactions among the 11 cell types under control and PE conditions, clearly showed a reduction in both metrics under PE conditions (Figure 7A and B). This suggests that signal transmission and interactions between cells may have been suppressed in PE, potentially affecting cellular function and overall tissue coordination. In particular, under PE conditions, endothelial cells, EVT, granulocytes, monocytes, myelocytes, and SCT exhibited reduced communication. These cell types play crucial roles in maintaining vascular stability, immune defense, and placental function, and their impaired communication can be linked to functional disruptions, thus playing a significant role in the pathogenesis of PE. However, the increased communication number and weight observed in CMP and VCT might reflect an activation or compensatory response under PE conditions. This increase in communication could have been a response to the reduction in other cellular interactions, aiming to maintain normal placental and hematological functions. Additionally, the bubble plot under PE conditions illustrated ligand-receptor interactions between these cell types, indicating a level of signaling and interaction and suggesting that these receptors and ligands play key roles in regulating cell behavior, maintaining tissue functions, and responding to environmental changes (Figure 7C).

Figure 7 Cellular communication analysis. (A) Networks of interactions between cell types (times). (B) Network of interactions between cell types (strength). (C) PE and control cell pairing with ligand receptors of other cells.

Identifying VCT as the Key Cell

VCT, EVT, and SCT are crucial components of the placenta. Dysfunction of these cells leads to abnormal development and malfunction of the placenta, which are key factors in the development of PE(12). Consequently, this study focused on the expression of biomarkers within these trophoblast layers, selecting cells exhibiting significant differences in biomarker expression as the key cells for further analysis. The results indicated that while only MIF and TLR4 expression differences were observed in EVT (P < 0.05), and MIF was the only biomarker differentiated in SCT (P < 0.05), significant differences in the expression of CD44, MIF, and PIK3R1 were found in VCT, all of which were also highly expressed under PE conditions (P < 0.05) (Figure 8A). Normally, VCT can differentiate into SCT; however, this process may have been disrupted during PE, thereby affecting placental function and maternal health. Thus, VCT was chosen as the key cell type for the subsequent pseudo-time analysis.

Figure 8 Identification of biomarkers. (A) The violin map of CD44, MIF, PIK3R1, TLR4 expression in differential cells. “NS” represented not significant, “*” represented p < 0.05, and “***” represented p < 0.001. (B) Distribution of biomarkers. (C) Distribution of biomarkers according to subgroups: red color represents the distribution of biomarkers, the left panel is the control group and the right panel is the PE group.

Additionally, the distribution of biomarkers across all 11 cell types was visualized, revealing that CD44 and TLR4 were predominantly found in granulocytes, whereas PIK3R1 was mostly present in macrophages (Figure 8B). Interestingly, an increase in MIF expression was observed in the VCT, EVT, and SCT under PE conditions (Figure 8C). As MIF regulates inflammatory and immune responses, its increased expression may reflect the placental cells’ reaction to abnormal pregnancy conditions. The upregulation of MIF may represent a biological response aimed at combating placental dysfunction and systemic inflammation caused by PE. In the future, modulation of MIF activity may offer a pathway to enhance placental function and improve pregnancy outcomes.

Differential Developmental Dynamics of VCT

Pseudo-time analysis was performed to determine the developmental trajectory of VCT. For this purpose, the VCT was re-clustered, identifying 13 distinct groups and 5 different stages of differentiation (Figure 9A). During the early (stage 1) and middle (stages 2 and 3) stages of VCT development, the number of cells was higher, whereas in the later stages (stages 4 and 5), the number of cells significantly decreased (Figure 9B and C). This decline could be attributed to the shift in VCT function from active proliferation and differentiation in the early and middle stages, which laid the foundation for placental maturation and development, to more specialized functions such as nutrient transfer and hormone synthesis in the later stages. A bifurcation in VCT development was noted as it progressed from the middle to later stages, differentiating into clusters 1, 5, and 9 (Figure 9D). This bifurcation typically indicates that a single-cell population begins to diverge into different cellular states or fates at this developmental point, highlighting the emergence of heterogeneity within the original cell population, which is substantial enough to steer cells along diverse trajectories.

Figure 9 Proposed time-series trajectory analysis. (A) Critical cell clustering results. (B) Difference in cell differentiation time (VCT), dark blue indicates an early stage of differentiation, while light blue indicates a later stage of differentiation. This can be used as a starting point for subsequent analysis. (C) Classes Occupied by Cell Differentiation (VCT), the different colors represent the taxa occupied by the cells. (D) Stages of Cell Differentiation (VCT), different colors represent different groups, 13 groups in total. (E) Cell differentiation in different samples. (F) Changes in the expression of biomarkers during differentiation of key cells.

Moreover, compared with the control, there was a notable decrease in the number of VCT in the early stages under PE conditions and an increase in the later stages, suggesting that PE might affect the normal growth and differentiation processes of placental cells (Figure 9E). The change in cell numbers could reflect the adaptive response of the placenta to an adverse pregnancy environment. The reduction in early stage VCT may have been due to insufficient oxygen and nutrient supply associated with PE, affecting the initial establishment and functional development of the placenta. Conversely, the increase in VCT numbers in the later stages could have been a compensatory response by the placenta to enhance its transport capabilities to support ongoing fetal development.

The expression levels of these biomarkers during VCT differentiation were monitored (Figure 9F). As VCT matured, CD44 showed an initial increase, followed by a decrease in expression; MIF exhibited a consistent upward trend, and PIK3R1 consistently decreased. No significant changes in TLR4 expression were observed.

Differential Expression of Biomarkers in PE and Controls

RT-qPCR results showed that the expression of CD44, TLR4, and MIF was upregulated in the PE group, while the expression of PIK3R1 was down-regulated (Figure 10).

Figure 10 The RT-qPCR results plots for D44, MIF, PIK3R1, and TLR4 genes. “ns” represented not significant, and “*” represented p < 0.05.

Discussion

PE is one of the main causes of maternal death and perinatal mortality; however, the specific mechanism remains unclear, and there is still a lack of effective prediction and intervention measures.⁶² Apoptosis is a form of programmed cell death and an important mechanism for maintaining tissue homeostasis. However, abnormal apoptosis of placental trophoblastic cells may lead to the occurrence of PE.¹¹ Through single-cell transcriptome technology to study apoptosis-related genes, in-depth exploration of potential biomarkers of PE revealed its potential mechanism of action, which has important clinical significance in improving the early diagnosis and effective management of PE.

In this study, all the data were sourced from public databases. Through the application of scRNA-seq analysis and differential expression analysis, coupled with five algorithms from the CytoHubba plugin, along with expression level analysis, four ARGs (CD44, MIF, PIK3R1, and TLR4) were identified as biomarkers for PE at the single-cell level. They were mainly associated with immunomodulatory-related pathways, with CD44 and TLR4 downregulation, and MIF and PIK3R1 showed increased expression in PE.

CD44 is a cell-surface glycoprotein that plays an important role in various biological processes, including inflammatory responses, tumorigenesis, and wound healing. It is considered one of the markers of cancer stem cells and is closely related to the aggressiveness and ability of tumors to metastasize, what’s more CD44 leads to the development of colorectal cancer by inducing apoptosis.^63,64 Placental trophoblasts are similar to tumor cells in that they have strong proliferative, migratory, and invasive abilities. Study has shown that increased CD44 positive expression was observed in maternal decidua cells and fibroblast cells close to the root villi.⁶⁵ Another study showed that the FKBPL-CD44 pathway appears to play a central role in the pathogenesis of PE and has potential for early diagnostic and therapeutic purposes in PE.⁶⁶ Taken together with previous studies, we hypothesized that CD44 may also play a role in PE apoptosis.

Macrophage Migration Inhibitory Factor (MIF) is a widely expressed cytokine that has pro-inflammatory effects, inhibits macrophage migration, and plays an important role in regulating immune responses, inflammatory responses and autoimmune disease.⁶⁷ In recent years, MIF has also been found to play a key role in a variety of pathological states such as tumorigenesis, metabolic diseases, and cardiovascular diseases.⁶⁸ It has been found that MIF is involved in the apoptotic pathway, and in a lipid metabolism disorder model, inhibition of the MIF pathway can reduce apoptosis of pancreatic islet β-cells and alleviate diabetes mellitus, whereas MIF inhibitors can also have a very good inhibitory effect on apoptosis in nasopharyngeal carcinoma.^69,70 A previous study showed that MIF levels were significantly increased in the maternal peripheral blood of the PE group.⁷¹ Similarly, we found that the expression of MIF in PE placental tissue samples was significantly higher than that in normal samples. It has been suggested that MIF can be used to regulate the inflammatory response of macrophages in the immune environment of PE and influence apoptosis.

The protein encoded by PIK3R1 is a regulatory subunits of phosphatidylinositol 3-kinase (PI3K), that regulates cell growth and differentiation.⁷² One study showed that, as an aging-related gene in PE, PIK3R1 may be related to different immune characteristics.⁷³ PIK3R1 is involved in the regulation of endometrial cell apoptosis.⁷⁴ Similar to our findings, another study revealed that PIK3R1 is the hub network and a key gene of PE, and needs more attention in future PE studies.⁷⁵

Toll-like Receptor 4 (TLR4) is a member of the Toll-like receptor family and is mainly expressed in immune cells, such as macrophages and dendritic cells. TLR4 activates the immune response by recognizing foreign pathogens and endogenous signaling molecules.⁷⁶ It plays a key role in the recognition of bacterial infections, the initiation of inflammatory responses, and immune regulation. Abnormal activation of TLR4 is associated with the development of multiple inflammatory, autoimmune, and chronic degenerative diseases.⁷⁷ TLR4 is directly involved in the regulation of apoptosis.⁷⁸ A previous study showed that, in early onset PE, the expression of TLR4 in the placental tissue was significantly higher, and the severity of PE was correlated with the degree of damage to the placental villi,⁷⁹ which is similar to our qPCR results.

All the above studies have reflected the key role of these four biomarkers in PE apoptosis. Subsequently, the present study integrated biomarkers to establish a nomogram for predicting PE, which showed good predictive ability and was superior to that of a single biomarker. This demonstrates the great potential of these biomarkers in clinical decision-making and is expected to provide new opportunities for the clinical diagnosis and treatment of PE.

To further explore the potential mechanisms of biomarkers in the development of PE, this study used GSEA enrichment analysis and found that four biomarkers were co-enriched in epithelial mesenchymal transition (EMT). EMT is a biological process in which epithelial cells lose their intercellular adhesion and polarity and gain the ability to migrate, invade, and transform into mesenchymal cells.⁸⁰ This process is critical for embryonic development, tissue remodeling, and cancer metastasis.⁸¹ In PE, EMT may be involved in placental formation and remodeling of the cell layer, especially in the invasion and migration of extra-villous trophoblastic cells (EVTs), which are essential for anchoring and functioning of the normal placenta. Abnormal EMT activity may lead to placental insufficiency and is associated with PE development of PE.⁸² Therefore, inhibition of EMT can reduce the abnormal activity of PE cells. However, the regulation of EMT by these biomarkers in PE remains unclear and an in-depth understanding of their regulatory mechanisms is of great significance for the treatment of PE.

Breakdown of maternal-fetal immune tolerance is one of the key etiologic factors in PE.⁸³ The present study also found that there were two different types of immune cells between the two groups: monocytes and resting NK cells. In fetal circulation, dominant monocyte intermediate subsets and increased inflammatory subsets are observed in PE compared to those in normal pregnancy.^84,85 In addition, NK Cell activation has been characterized and is known to play major causative roles in the maternal inflammatory response in PE.^86,87 In addition, the activity of APC co-stimulation, Checkpoint, T cell co-stimulation, and MHC class I pathways was reduced in the PE group, indicating immunosuppression in patients with PE. Immune checkpoints are molecules in the immune system that control and regulate the activity of immune cells by inhibiting signaling and regulating the immune response, thereby helping to maintain the immune system to attack pathogens while avoiding damage to normal tissue.⁸⁸ Because of immune dysfunction in PE, immune checkpoint inhibitors are emerging as new therapeutic tools for PE, where CD44 and T cell co-stimulation can effectively improve the immune efficiency of vaccines.^89,90 MHC I and TLR4 play important roles in endogenous antigen presentation and recognition of exogenous pathogens. Both synergistically participate in the immune response while avoiding excessive damage to the host tissues.⁹¹ Based on the differential immune cells and different immune pathways in PE, it has been suggested that targeted biomarkers can be used to regulate the immune infiltration of PE and that T cell co-stimulation or MHC class I expression can be regulated by biomarkers to enhance the immune response of PE. This study aimed to provide new immunotherapy for the clinical treatment of PE and to improve the outcomes of PE.

TLR4 predicts drug isoflavones. Isoflavones are natural plant compounds whose biological activity is mainly related to estrogen receptors and plays an important biological role in cancer, cardiovascular disease, oxidative stress, and other diseases through TRL4.⁹² Puerarin (Pue) is an isoflavone that is derived from puerarin. Pue can improve trophoblast movement by inhibiting the PE-related inflammatory response, and is speculated to have a potential therapeutic effect on PE.⁹³ This study established a drug-biomarker-immune cell regulatory network that provides a broader vision for PE immunotherapy.

PE is a complex condition that affects pregnancy and involves multiple cell types including VCT, these includeVCT SCT, and EVT. These cell types play key roles in the formation and function of placenta.^94–96 Single-cell analyses confirmed the importance of these three cell types. Finally, in the RT-qPCR experiments, all three biorepresentations showed the opposite trend, except for MIF, which showed the same expression trend as that in the dataset. Combined with single-cell analysis, this study also found that biomarkers showed dynamic expression in VCT cells. Taken together with previous studies, we speculate that this may be due to the period in which the current PE samples were collected, as well as the number of samples, and that more samples will need to be collected in the future for in-depth validation and study of biomarkers.

Although the nomogram constructed in this study shows great potential in predicting preeclampsia (PE), further in-depth exploration is needed on how to integrate it into the existing clinical workflow. Currently, there are many limitations in the early prediction methods for PE in clinical practice, such as poor accuracy and insufficient specificity,⁹⁷ which makes the early and accurate diagnosis of PE and timely intervention challenging. From a practical application perspective, integrating the nomogram into the clinical workflow first requires addressing the issue of data acquisition. The use of the nomogram relies on accurate and comprehensive clinical data, including patients’ basic information, medical history, and laboratory test indicators. Clinicians need to ensure that they can efficiently and accurately collect these data during daily diagnosis and treatment and input them into the nomogram model.⁹⁸ For example, the electronic medical record system can be optimized to automatically extract and integrate relevant data and directly connect to the nomogram prediction model, reducing errors and time costs caused by manual operations.⁹⁹ In terms of overcoming the shortcomings of existing early prediction methods, the nomogram has unique advantages. It comprehensively considers the expression levels of multiple biomarkers related to PE (such as CD44, MIF, PIK3R1, and TLR4), providing more comprehensive and accurate predictive information compared to the detection of a single indicator.¹⁰⁰ However, in practical application, the reliability and stability of the nomogram in different populations and clinical environments still need to be further verified. Multi-center, large-sample clinical studies can be conducted to collect data from patients in different regions and races, evaluate the predictive efficacy of the nomogram, and optimize and adjust the model based on the results.¹⁰¹ In addition, the acceptance of the nomogram by doctors and patients also needs to be considered in clinical application. Doctors should be trained and patients should be well informed.¹⁰² In conclusion, although the nomogram provides a new tool for the early prediction of PE, its wide application in clinical practice still requires in-depth research and practical exploration in aspects such as data acquisition, model validation, doctor training, and patient education.

Through single-cell transcriptome sequencing and bioinformatics analysis, four ARGs (CD44, MIF, PIK3R1, and TLR4) were identified as biomarkers, and the same expression trend was observed in the dataset and PE clinical samples. These biomarkers are involved in immune infiltration, oxidative phosphorylation, and other processes in PE and are expected to become potential markers of PE. However, clinical applications require more data support from samples. The effects of these key genes on PE were further studied in animal and cell experiments.

Conclusion

PE is a hypertensive multisystem disorder that causes significant fetal-maternal mortality and morbidity worldwide. Apoptosis occurs in PE and is associated with immune cells and pathways. According to the present study, four ARGs, CD44, MIF, PIK3R1, and TLR4, offer new insights into PE’s cellular mechanisms of PE and may serve as potential biomarkers for PE. This study opens new avenues for future research on PE’s pathogenesis of PE and treatment options. Nevertheless, this study has some limitations. The sample size used in the current research is relatively small, making it difficult to cover the genetic, physiological, and environmental differences among individuals, which limits the generalizability of the research conclusions and affects the accuracy and stability of the nomogram model predictions. Future research will be conducted on a large scale, multi-center, and will include different patient populations; in-depth functional studies will be carried out using in vitro cell experiments and in vivo animal models; the interactions between these genes and other related factors of PE will be explored in depth to gain a more comprehensive understanding of the disease mechanism, aiding in the development of more effective treatment strategies.