The GSE4707 and GSE10588 datasets were downloaded from the GEO database (http://www.ncbi.nlm.nih.gov/geo/). GSE4707 was based on GPL1708, with an Agilent-012391 Whole Human Genome Oligo Microarray G4112A (Feature Number version), and placental biopsies were obtained from ten patients with preeclampsia and four women with normal pregnancies during Caesarean section. Chorionic tissue was dissected from a standardized location approximately 2 cm beside the umbilical cord insertion, from the middle layer of the placenta midway between maternal and fetal surfaces. In the GSE4707 dataset, there was no statistical difference in gestational age and weight in patients with preeclampsia compared with normal mothers, but neonatal weight in the preeclampsia group was significantly lower than neonatal weight in the normal group. GSE10588 was based on GPL2986 using ABI Human Genome Survey Microarray Version 2, and placental biopsies were obtained from 17 patients with preeclampsia and 26 women with normal pregnancies. A central area of chorionic tissue was dissected, and the maternal deciduas and amnionic membranes were removed. Scientists then dissected 1-cm-thick sections of placental villi from the central area between basal and chorionic plates. In the GSE10588 dataset, there were no statistical differences in maternal age, BMI, but patients with preeclampsia had lower gestational age and higher cesarean delivery rates compared to normal mothers. The probe in each data file was modified to a gene symbol according to its annotations in the probe file. If the same gene symbol corresponded to various probes, the diameter of the probe was assessed as the final expression value. Then, the results were merged for further analysis. Simultaneously, the function of the “SVA” package of R version 4.2.0 was used to eliminate the batch effect (Leek et al., 2012). We used the GSE160888 dataset according to Agilent-045997 Arraystar human lncRNA microarray V3, which included placental samples from four preeclampsia and four control cases. The placenta specimen was resected from the middle of the villous lobule. We also used GSE96985 based on Agilent-078298 human ceRNA array V1.0 4 × 180K, which included placental samples from four preeclampsia and three control cases. Using the method mentioned above, both datasets were combined as a validation cohort. The GSE48424 dataset, which was based on Agilent-014850 Whole Human Genome Microarray 4 × 44K G4112F and includes whole blood samples from 18 preeclampsia and 18 control cases, was also downloaded.

The two datasets were combined into a metadata queue, and the SVA package was used to eliminate the batch effect of the two datasets. The R limma package (http://www.bioconductor.org/) was used for background correction, normalization, and differential expression between array analyses. We used a false discovery rate-adjusted sample (false discovery rate adjusted; p < 0.05) and | log2 Fold Change | > 1 as the threshold point for DEGs.

To determine the main biological properties of DEGs, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed. The R packages “ggplot2” (Wickham, 2011), “enrich Plot”, and “clusterProfiler” (Yu et al., 2012) were used to create GO and KEGG enrichment plots. Disease ontology (DO) enrichment analysis was also conducted on DEGs using “clusterProfiler” and the “DOSE” software package in R (Yu et al., 2012; Yu et al., 2015). Through Gene Set Enrichment Analysis (GSEA), the most notable functional differences between the preeclampsia and control groups were verified (Subramanian et al., 2005). The study set c2.cp.go.v7.0. symbols, and the GMT from the Molecular Marker Database (MSigDB) were used as reference gene sets (Liberzon et al., 2015).

We used three machine-learning algorithms to screen for important biomarkers of preeclampsia. The least absolute shrinkage and selection operator (LASSO) is a regression analysis algorithm that uses regularization to increase prediction accuracy. In the “glmnet” software package in R, the LASSO regression algorithm was used to select preeclampsia genes that were notably related to normal samples (Engebretsen and Bohlin, 2019). Support vector machine (SVM) is a popular machine-learning technique for classification and regression (Nedaie and Najafi, 2018). We used “E1071,” “Kernlab,” and “caret” to build the SVM model. To avoid overfitting, we determined the optimal genes from the metadata queue using the Recursive feature elimination (RFE) algorithm. The LASSO and SVM were used to select overlapping genes, and the level of candidate genes was verified on the GSE160888 and GSE96985 datasets. In addition, we established a random forest (RF) model in the “RandomForest” package in R as a training model to forecast the occurrence of preeclampsia. Moreover, ntrees and mtry were set at 100 and 3, respectively. Next, a rosette model was built with the “RMS” package in R to forecast the prevalence of patients with preeclampsia. Calibration curves were used to determine the agreement between predicted and actual values. A clinical influence curve was drawn using decision curve analysis (DCA) to evaluate whether model-based decision-making is beneficial to patients (Iasonos et al., 2008). We used the “limma”, “ggplot2”, “ggpubr”, and “ggExtra” packages in R to calculate and plot the correlations between the four characteristic genes.

To validate the predictive value of previously screened biomarkers, receiver operating characteristic (ROC) curves were built via the “pROC” package based on mRNA expression data using 30 preeclampsia and 27 control samples. The area under the ROC curve (AUC) was used to judge the diagnostic efficiency of preeclampsia and control samples, which was further verified using the GSE160888, GSE96985, and GSE48424 data files.

The CIBERSORT bioinformatics (https://cibersortx.stanford.edu/) algorithm was used to analyze immune cell infiltration from the preeclampsia gene expression profiles. Putative immune cell abundance was predicted using a reference set of 22 immune cell subtypes with 1000 permutations (LM22; Newman et al., 2015). Correlation analysis and visualization of 21 infiltrating immune cells were performed using the R package “Corrplot.” The Vioplot software package was used to draw a violin diagram to observe the difference in immune cell infiltration between patients with preeclampsia and normal pregnant women.

Pearson’s correlation analysis was performed using R to determine the correlation between the identified gene biomarkers and the level of infiltrating immune cells. The resulting associations were visualized using the graphical techniques of the “ggplot2” package.

Twenty-one pregnant women were selected from the International Peace Maternal and Child Health Hospital affiliated to Shanghai Jiao Tong University School of Medicine. Nine cases were diagnosed as preeclampsia and twelve cases were normal pregnant women. Preeclampsia was diagnosed according to the ACOG Practice Bulletin (The American College of Obstetricians and Gynecologists, 2020). All participants were singleton pregnant women without other diseases affecting blood pressure, such as hyperthyroidism, Cushing’s syndrome, and pancreatitis; serious dysfunction of the heart, liver, and kidney; and acute complications, such as diabetic ketoacidosis. All pregnant women provided signed informed consent, and this study was approved by the Ethics Committee of International Peace Maternal and Child Health Hospital Affiliated to Shanghai Jiao Tong University School of Medicine [Approval No.: GKLW-2017-81].

General information on the pregnant women, including age at delivery and family history of hypertension, was collected through face-to-face interviews. The pre-pregnancy body-mass index was calculated from the height measured by nurses and the pre-pregnancy weight reported by the pregnant women. Blood pressure was measured during the second trimester.

The amniotic membrane was removed within 5 min after delivery of the placenta. We dissected the middle layer of the placenta from the middle of the maternal and fetal surfaces, followed by washing with enzyme-free water to remove blood. Then, the water was removed using filter paper, and the placenta was quickly placed in an external rotating freezing tube, flash-frozen in liquid nitrogen, and stored at −80°C for later use.

The total RNA of the placenta was extracted using RNAiso Plus reagent (9109, Takara, Shiga, Japan) according to the manufacturer’s instructions. cDNA was synthesized from RNA using the RT Reagent Kit and gDNA Eraser (RR047A, Takara, Shiga, Japan). qPCR was performed on the QuantStudio 7 Flex system (Life Technologies, Carlsbad, CA, United States). Three replicates of each sample were analyzed. To quantify the relative mRNA expression, data were normalized to the expression level of β-actin. Primer sequences are shown in Table 1.

R and SPSS26.0 were used for statistical analyses. Student’s t-test and Mann-Whitney test were used for continuous, normally and non-normally distributed data, respectively. The “glmnet” and E1071 packages in R and ROC curve analysis were used for the LASSO regression algorithm and the diagnostic efficacy of the selected biomarker analysis, respectively. Pearson’s correlation coefficient was used to study the relationship between gene biomarker expression and immune cell infiltration. For quantitative data, distribution was described in terms of mean ± standard deviation, and an independent sample t-test was used for comparison between groups. The median M [P25, P75] was used to describe data that did not conform to normal distribution, and the rank sum test of independent samples was used for comparison between groups. Enumeration data were described in terms of the number of cases (%), and the Fisher’s exact probability method was used for comparison of differences between groups. All statistical analyses were bilateral, and p < 0.05 was considered to indicate significance.

Two GEO datasets (GSE4707 and GSE10588) were downloaded, which together included 27 patients with preeclampsia and 30 normal pregnant women. After eliminating batch sub-effects, the “limma” package was used to determine the DEGs. We created heat maps with the screened differential genes, with red representing upregulated genes in PE patients and blue representing downregulated genes (Figure 2A). We also plotted volcano plots, with red representing upregulated genes in PCOS patients and green representing downregulated genes (Figure 2B). Seventy-three DEGs were identified, of which 56 were significantly upregulated and 17 significantly downregulated.

GO analyses showed that DEGs were predominantly enriched in reproductive structure development, reproductive system development, and hormone metabolic processes (Figure 3A). KEGG analyses indicated that DEGs were predominantly enriched in cytokine-cytokine receptor interaction, cell adhesion molecules, and the interleukin (IL)-17 signaling pathway (Figure 3B). DO pathway enrichment analyses suggested that diseases enriched by DEGs were largely associated with endocrine system diseases, preeclampsia, and reproductive system diseases (Figure 3C). GSEA results of preeclampsia showed that the enriched pathways dominated myeloid cell homeostasis and cadherin binding (Figure 3D).

Two different algorithms were used to identify promising biomarkers of preeclampsia. First, 12 genes were identified as diagnostic biomarkers of preeclampsia using the LASSO algorithm based on the DEGs (Figure 4A). Then, the SVM-recursive feature elimination algorithm was used to screen ten characteristic genes from the DEGs (Figure 4B). Four overlapping features (leptin [LEP], SAM and SH3 domain containing 1 [SASH1], RAB6C, and fms-like tyrosine kinase receptor-1 [FLT1]) were selected using the two algorithms (Figure 4C). We also constructed an RF tree and screened out two genes with a score greater than two according to their importance score, namely LEP and SASH1 (Figures 4D, E). We then used the GSE160888 and GSE96985 datasets to analyze their levels. FLT1 levels in the placenta of patients with preeclampsia were significantly upregulated (p = 0.014; Figure 5A). Furthermore, the expression of LEP in the tissues of patients with preeclampsia was higher than that in the control group and that of RAB6C was lower, but no significant differences were observed (Figures 5B, C). Furthermore, the SASH1 levels in the placenta of patients with preeclampsia were significantly upregulated (p = 0.028; Figure 5D). We then used the GSE48424 dataset, which contained 18 blood samples from patients with preeclampsia and 18 blood samples from control patients, to analyze the expression of these four genes. FLT1 expression in the blood samples of patients with preeclampsia was higher than that in control patients, but the difference was not significant (Figure 5E). The expression of LEP in the blood samples of patients with preeclampsia was significantly higher than that in control patients (p = 0.047; Figure 5F); the expression of RAB6C was lower in the blood of patients with preeclampsia (Figure 5G). This was consistent with the results validated in the placenta. However, SASH1 showed the opposite results in the placenta and blood. The expression of SASH1 was significantly increased in the placenta of patients with preeclampsia but decreased in the blood (Figures 5D, H).

We examined the expression of four genes in nine patients with preeclampsia and 12 normal pregnant women. Basic information on the patients is shown in Table 2. The expression of SASH1 was significantly increased in patients with preeclampsia both in the predicted and qPCR experiment results (Figures 5D, L). However, we predicted that FLT1 expression would be significantly increased in preeclampsia (Figure 5A), but there was no significant difference in the qPCR results (Figure 5I). Both the predicted results and the experimental verification showed that the expression of LEP in patients with preeclampsia was higher and that of RAB6C was lower than those in normal pregnant women, but there was no significant difference (Figures 5B, J).

We analyzed the correlations among the expression of FLT1, LEP, SASH1, and RAB6C in GSE4707 and GSE10588. LEP was positively correlated with FLT1 (R = 0.87, p < 2.2E-16; Figure 5M), RAB6C (R = 0.65, p = 4.9E-08; Figure 5N), and SASH1 (R = 0.73, p = 1.6E-10; Figure 5O). Furthermore, SASH1 was positively correlated with FLT1 (R = 0.64, p = 6.6E-08; Figure 5P) and RAB6C (R = 0.57, p = 5.2E-06; Figure 5Q). Finally, RAB6C was positively correlated with FLT1 (R = 0.64, p = 8.3e-08; Figure 5R).

A nomogram model for predicting the prevalence of preeclampsia based on the two candidate genes was constructed using the ‘RMS’ package (Figure 6A). The calibration curve showed that the prediction ability of the nomogram model was optimized (Figure 6B). The red line in the DCA curve remained above the gray and black lines between 0 and 1, indicating that decisions based on the nomogram model may be beneficial to patients with preeclampsia (Figure 6C). The clinical influence curve showed that the nomogram model had a good predictive ability (Figure 6D).

The four biomarkers showed good diagnostic value for preeclampsia and control samples in the GSE4707 and GSE10588 datasets (Figure 7A). The AUC of FLT1 was 0.922 (95% confidence interval [CI] 0.839–0.986) and that of LEP was 0.967 (95% CI 0.915–0.999). The AUC of SASH1 was 0.954 (95% CI 0.890–1.000) and that of RAB6C was 0.937 (95% CI 0.870–0.984). The GSE160888 and GSE96985 datasets verified its recognition ability, with an AUC of 0.875 (95% CI 0.571–1.000) for FLT1 and 0.714 (95% CI 0.411–0.946) for LEP. The AUC of RAB6C was 0.589 (95% CI 0.286–0.875) and that of SASH1 was 0.848 (95% CI 0.607–1.000), indicating that the characteristic biomarkers had diagnostic ability (Figure 7B). The GSE48424 dataset showed verified recognition ability, with an AUC of 0.515 (95% CI 0.327–0.728) for FLT1 and 0.694 (95% CI 0.494–0.858) for LEP. The AUC of RAB6C was 0.688 (95% CI 0.497–0.855), and that of SASH1 was 0.577 (95% CI 0.370–0.759), indicating that the characteristic biomarkers had diagnostic ability (Figure 7C).

Only activated NK cells and resting dendritic cells were statistically significant (p < 0.05; Figure 8A), whereas the other cell types did not differ significantly between the placenta of normal pregnant women and patients with preeclampsia. The correlation between infiltrated immune cells is displayed in Figure 8B; Table 3.

FLT1 was positively correlated with the abundance of Tregs (r = 0.28, p = 0.034) but negatively correlated with that of activated NK cells (r = −0.41, p = 0.001; Figure 9A). RAB6C was negatively correlated with the abundance of M1 macrophages (r = −0.30, p = 0.025), resting dendritic cells (r = −0.32, p = 0.014), and activated NK cells (r = −0.43, p < 0.001; Figure 9B). LEP was negatively correlated with the abundance of activated NK cells (r = −0.32, p = 0.016; Figure 9C). SASH1 was negatively correlated with the abundance of activated NK cells (r = −0.28, p = 0.034) and resting dendritic cells (r = −0.31, p = 0.018; Figure 9D).