First, we obtained the matrix files of GSE24129, GSE54618, GSE60438, and GSE75010 from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/gds) (Table 1). On the basis of their probe annotation files, the probes were mapped to gene symbols in each data set. For multiple probes corresponding to the same gene, the average expression value of the gene was calculated to represent the gene expression level. After merging the three data sets into a metadata cohort and removing the batch effect, the “SVA” package of R software was applied (Leek et al., 2012). The background correction and normalization of raw data were processed by the limma package of R (http://www.bioconductor.org/), where genes with |log fold change (FC)|>0.5 and adjusted p < 0.05 were defined as DEGs. GSE75010 was used to determine and validate functions of significant DEGs (Table 1).

A total of 15 differentially expressed genes (DEGs) were analyzed using the R language with the clusterProfiler, org.Hs.eg.db, enrichplot, and ggplot2 packages via Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses. p < 0.05 was considered to denote significant enrichment. The R package “clusterProfiler” and “DOSE” were applied to conduct the Disease Ontology (DO) enrichment analyses on DEGs. The Gene Set Enrichment Analysis (GSEA) was used to analyze the association between the risk score and hallmarks. The GSEA was applied to recognize the most important feature among the PE and normotensive pregnancy groups. The “c5.go.v7.4.symbols.gmt” was applied as the reference gene set from the Molecular Signatures Database (MSigDB). The gene set was considered to be significant at a p < 0.05.

We predicted disease status via two machine-learning algorithms. To increase the prediction accuracy, the least absolute shrinkage and selection operator (LASSO)—a regression-based analysis method that scrutinizes variable selection and regularization in PE models—was used. The R package “glmnet” was applied to conduct LASSO regression analysis on identification of valuable DEGs related to the discrimination of PE and normotensive control. The support vector machine (SVM) is an efficient and widely applied supervised machine-learning algorithm for disease classification and regression tasks (CORTES and Vapnik, 1995). Hence, conjugated LASSO and SVM-RFE were used to screen the overlapping genes, which were further verified in the GSE75010 data set.

The mRNA expression data were obtained from 73 PE and 85 normotensive pregnancy samples from the GEO database, which were used to generate the receiver operating characteristic (ROC) curves to determine the predictive value of biomarkers. The area under the ROC curve (AUC) was used to ascertain the diagnostic capability of biomarkers in distinguishing PE earlier from normotensive samples and was also verified by the GSE75010 data set.

To build a PPI network, we employed “LEPTIN” in the “protein name” module and “Homo sapiens” in the organism module to search from the STRING website (https://string-db.org/). We set the key parameters as follows: meaning of network edges (“confidence”), the minimum required interaction score [“highest confidence (0.900)”], and the maximum number of interactors to show (“no more than 10 interactors”) in the first shell. The STRING tool was used to perform KEGG and GO molecular function analyses of LEP-related genes.

From the gene expression profiles in PE, the CIBERSORT (http://cibersort.stanford.edu/) algorithm was used to quantify the proportion of immune infiltration cells. We downloaded a gene signature matrix with interpretation, known as the 22 kinds of immune cells (LM22) with 1,000 permutations, from the webpage of CIBERSORT to estimate the putative abundance of immune cells (Newman et al., 2015). The R package “corrplot” was applied to conduct the correlation and visualization of LM22. The R package “vioplot” was used to visualize the differences in immune cells between the PE and normotensive control groups. Pearson’s correlation analysis was used to explore the screened diagnostic biomarker relation to the levels of immune infiltration cells. The chart technique with R package “ggplot2” was applied to visualize the aforementioned result.

A total of 32 paraffin-embedded PE and 41 normotensive specimens were diagnosed at the Second Affiliated Hospital of Fujian Medical University (Fujian, China) from August 2017 to September 2021. The research was approved by the Research Ethics Committee of the Second Affiliated Hospital of Fujian Medical University prior to the study.

IHC staining was performed as previously described (Chen et al., 2016). The primary antibody was anti-leptin (Servicebio, Wuhan). Leptin staining intensity ratios were scored as follows: negative = 0, light yellow = 1, brownish yellow = 2, or tan = 3. The staining cells were scored as follows: less than 1/3 = 1, between 1/3 and 2/3 = 2, or more than 2/3 = 3. The final score for leptin expression was calculated by multiplying the two scores. Slides were divided into low- and high-expression groups, corresponding to scores of <3 or ≥3, respectively. The histopathological diagnosis of the patients was established by two pathologists specialized in obstetrics and gynecology.

Total RNA was extracted from placental tissue immediately after normal labor and cesarean section utilizing TRIzol (TaKaRa, Japan), and then, cDNA was prepared according to the protocol (TaKaRa, Japan). The detailed procedure is presented in the Supplementary Methods. GAPDH was used as an internal reference and the relative mRNA expression of leptin was calculated by the 2^−ΔΔCT method. qRT-PCR for each sample was repeated in three independent experiments. The primer sequences are shown below:

GAPDH

Forward: 5′-CAT​GTT​CGT​CAT​GGG​TGT​GAA​CCA-3′,

LEPTIN

Forward: 5′-AAC​GTG​ATC​CAA​ATA​TCC​AAC​G-3′,

Total protein was extracted from each 50 mg placenta sample using RIPA lysis buffer (CW Biotechnology, Beijing China) with protease (Solarbio, China), DNA enzyme inhibitor (Solarbio, China), and phosphatase inhibitors (CW Biotechnology, Beijing China). Western blotting was performed as described previously (Guo et al., 2016). The following antibodies were utilized: anti-leptin (Servicebio, China) and anti-GAPDH (Cell Signaling Technology, United States). The secondary antibody was anti-rabbit IgG (Cell Signaling Technology, United States). Gray values were measured with ImageJ.

We used the R software (v.4.1.1) for all statistical analyses. The Mann–Whitney U test was used to compare the PE and normotensive control groups. The LASSO regression analysis, SVM algorithm, ROC curve analysis, Pearson’s correlation, and unpaired t-test were used as described above. p < 0.05 was considered to indicate statistically significant differences for all analyses.

The analysis procedure used in the study is shown in Figure 1. We downloaded the transcriptome RNA-seq data from the GEO database. We identified DEGs between PE and normotensive control groups. We conducted GO, KEGG, DO, and GSEA analyses on DEGs. Conjugated LASSO and SVM-RFE were used to screen the overlapping candidate genes; the ROC curve was applied to determine the predictive value of the biomarkers, which were further validated in the GSE75010 data set. The CIBERSORT algorithm was used to calculate the compositional patterns of LM22 in PE. Correlation analysis was performed among the diagnostic markers and infiltrating immune cells.

This study involved three data sets from the GEO database (GSE24129, GSE54618, and GSE60438) and included a total of 73 PE and 85 normotensive control pregnancy samples. We obtained a total of 15 DEGs when comparing PE and normotensive controls (Figures 2A, B). Of these genes, four were remarkably downregulated and 11 were markedly upregulated (Figure 2B).

The results from the GO analysis showed that the DEGs were closely related to hormone-related GO terms, such as hormone activity, steroid metabolic process, and positive regulation of hormone secretion, and to receptor ligand–related GO terms, particularly receptor ligand activity, signaling receptor activator activity, and peptide hormone receptor binding (P < 0.05, Figure 3A). In addition, KEGG analysis showed enrichment of the GnRH secretion, neuroactive ligand–receptor interaction, cAMP signaling pathway, and cytokine–cytokine receptor interaction (P < 0.05, Figure 3B). The analysis of DO enrichment showed that DEGs were mostly related to cell-type benign neoplasm, infertility, PE, and kidney failure (Figure 3C). The GSEA results revealed that the enriched functions mainly enrolled activation of immune response, adaptive immune response, and cell chemotaxis in control (Figure 3D, Supplementary Table S1); and glycosylation, multicellular organism process, and hormone activity in PE (Figure 3E, Supplementary Table S2). These results conclusively show that immune response and hormone secretion have vital roles in PE. A PPI network with 15 DEGs was obtained. The STRING tool identified five binding proteins. LEP, LHB, PDK4, SP1, and KRT19 are shown in Figure 3F.

We applied the LASSO and SVM-RFE algorithm to selective potential biomarkers. The seven DEGs were identified as diagnostic biomarkers using LASSO regression for PE (Figure 4A). The five DEGs were verified by applying the SVM-RFE (Figure 4B). Between the two algorithms, the seven overlapping candidate genes (LEP, PDK4, SPP1, NAPSB, CGB5, LTB, and FCN1) were ultimately screened (Figure 4C). Additionally, to produce more reliable and accurate DEGs, verification of the seven DEGs’ expression levels was conducted using the GSE75010 data set. The CGB5, LEP, LTB, and NAPSB expression levels in PE samples were remarkably higher than in that of the normotensive group (Figures 5A, C–E; p < 0.05). The PDK4 and SPP1 expression levels in PE tissue were significantly lower than those in the normotensive group (Figures 5F, G; p < 0.05). However, the FCN1 expression showed no significant difference between the two groups (Figure 5B). Hence, we next explored the potential value of the diagnostic model combination of the six identified DEGs by applying a logistic regression algorithm.

Both LEP (AUC = 0.712) and PDK4 (AUC = 0.718) showed a good diagnostic value for the early diagnosis of PE (Figure 6A). Furthermore, a forceful discerning capacity was verified in the GSE75010 data set with an AUC of 0.850 in LEP (Figure 6B), showing that the LEP gene had a higher diagnostic capacity. We assessed the expression of LEP across PE and normotensive tissues via immunohistochemistry and found that high expression of LEP was associated with PE (Figure 6C; p < 0.05). For further clinical validation, we assessed the expression of LEP using qRT-PCR and revealed that high expression of LEP was associated with PE (Figure 6D; p < 0.05). We performed Western blotting to assess that LEP is upregulated in PE and found that high expression of LEP was associated with PE (Figure 6E; p < 0.05). The aforementioned results indicate that the LEP gene had a higher diagnostic capacity.

A PPI network with 14 nodes was obtained. The STRING database analysis identified 10 highest leptin-binding proteins. As shown in Figure 7, GCG, IAPP, LEPR, GHRL, STAT3, PPARG, SOCS3, NPY, JAK2, and PTPN1 were predicted to have the most powerful interactions with leptin. And then we performed the GO molecular function and KEGG analyses of the proteins predicted via the STRING tool. The GO_MF analysis showed significant differences in the relevance of the signaling receptor binding and hormone activity (Table 2). KEGG analysis showed significant differences in association with the adipocytokine and JAK-STAT signaling pathways (Table 3). The KEGG analysis revealed that STAT3, SOCS3, LEPR, JAK2, and NPY were associated with leptin in the adipocytokine signaling pathway, and STAT3, SOCS3, LEPR, and JAK2 were associated with leptin in the JAK-STAT signaling pathway, suggesting that leptin may interact with these proteins to activate both adipocytokine and JAK/STAT signaling in PE.

To further verify the relationship between the LEP gene and immune cell infiltration, we first applied the CIBERSORT algorithm to determine the proportions of the 22 types of infiltrating immune cells in the PE and control samples (Figures 8A, B; Supplementary Table S3). Next, the component of immune cells in PE vs. normotensive group was explored. The ratio of monocytes in PE was significantly lower than that in the normotensive control (P < 0.001). However, the ratio of CD4⁺ resting memory T cells (p = 0.015) in PE was significantly higher than that in normotensive controls (Figure 8C). Furthermore, we studied the relationship between the LEP gene and infiltrating immune cells. LEP was positively correlated with gamma delta T cells (r = 0.237, p < 0.05), M0 macrophages (r = 0.224, p < 0.05), memory B cells (r = 0.192, p < 0.05), and regulatory T cells (r = 0.181, p < 0.05), while being negatively correlated with resting CD4 memory T cells (r = −0.208, p < 0.05), M0 macrophages (r = −0.205, p < 0.05), and activated NK cells (r = −0.200, p < 0.05) (Figure 8D). The impact of the LEP gene was supported by these results regarding immune activity.