https://orcid.org/0009-0003-2059-5160
Figures
Abstract
Objective
Preeclampsia (PE) is a severe complication of unclear pathogenesis associated with pregnancy. This research aimed to elucidate the properties of immune cell infiltration and potential biomarkers of PE based on bioinformatics analysis.
Materials and methods
Two PE datasets were imported from the Gene ExpressioOmnibus (GEO) and screened to identify differentially expressed genes (DEGs). Significant module genes were identified by weighted gene co-expression network analysis (WGCNA). DEGs that interacted with key module genes (GLu-DEGs) were analyzed further by Kyoto Encyclopedia of Genes and Genomes (KEGG) and Gene Ontology (GO) analyses. The diagnostic value of the genes was assessed using receiver operating characteristic (ROC) curves and protein-protein interaction (PPI) networks were constructed using GeneMANIA, and GSVA analysis was performed using the MSigDB database. Immune cell infiltration was analyzed using the TISIDB database, and StarBase and Cytoscape were used to construct an RBP-mRNA network. The identified hub genes were validated in two independent datasets. For further confirmation, placental tissue from healthy pregnant women and women with PE were collected and analyzed using both RT-qPCR and immunohistochemistry.
Results
A total of seven GLu-DEGs were obtained and were found to be involved in pathways associated with the transport of sulfur compounds, PPAR signaling, and energy metabolism, shown by GO and KEGG analyses. GSVA indicated significant increases in adipocytokine signaling. Furthermore, single-sample Gene Set Enrichment Analysis (ssGSEA) indicated that the levels of activated B cells and T follicular helper cells were significantly increased in the PE group and were negatively correlated with GLu-DEGs, suggesting their potential importance.
Citation: Mao Y, Li X, Ren R, Yuan Y, Wang L, Zhang X (2024) Identification of hub glutamine metabolism-associated genes and immune characteristics in pre-eclampsia. PLoS ONE 19(5): e0303471. https://doi.org/10.1371/journal.pone.0303471
Editor: Jinhui Liu, The First Affiliated Hospital of Nanjing Medical University, CHINA
Received: December 19, 2023; Accepted: April 12, 2024; Published: May 8, 2024
Copyright: © 2024 Mao et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Data Availability: All relevant data are within the paper and its Supporting Information files.
Funding: This study was supported by grants from the Natural Science Foundation of Gansu Province(21JR7RA608).The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
Introduction
Preeclampsia (PE) is a progressive multisystem disease that is unique to pregnancy and significantly threatens the life and health of mothers and children. The average global incidence of PE is 3‒8% [1]. The PE-associated pathophysiological changes develop in early pregnancy; however, hypertension and proteinuria are observed after 20 weeks with or without progressive multi-organ damage. Currently, clinical therapy is limited to symptomatic treatment, and there are no effective preventive measures to improve the health of the mother and child. Emergency cesarean sections are often performed, resulting in serious iatrogenic preterm birth and complications associated with cesarean sections [2]. Therefore, it is important to identify markers of PE that will allow its early identification, thus reducing the maternal mortality rate and improving the neonatal survival rate.
Immunological factors are thought to be significantly associated with the etiology and pathogenesis of PE. Immune homeostasis is essential for the maintenance of pregnancy, and studies have indicated that abnormal activation of the immune system during pregnancy might be significantly associated with PE development [3], such as the activation of monocytes and neutrophils, abnormal cytokine production, and imbalances between regulatory T cells and helper T cells [4].
Glutamine (GLU) is the most abundant amino acid in the circulation and is rapidly utilized by tumor cells. It is frequently utilized as a citrate source for lipid production in reductive carboxylation or for maintaining TCA flux in cellular aerobic glycolysis. Furthermore, glutaminolysis assists the survival of proliferating cells by reducing oxidative stress and maintaining mitochondrial membrane integrity. It also contributes to the energy supplies of immune and tumor cells. However, inflammatory immune cells such as macrophages that act against the tumors do not utilize GLU metabolism. In comparison with naïve macrophages, M2 macrophages are dependent on GLU, whereas inhibition of GLU metabolism induces the pro-inflammatory M1 macrophage phenotype [5]. Thus, GLU metabolism is a possible target for the phenotypic conversion of tumor-associated macrophages (TAMs) from the M2 to the M1 phenotype, which increases the anti-tumor responses of immune-mediated inflammation. Additionally, GLU metabolism is associated with Th1 cell differentiation and effector T cell activation [6]. These data indicate that targeting of GLU metabolism can remodel the tumor microenvironment (TME) and enhance the efficacy of immunotherapy. It has been reported that suppression of GLU metabolism substantially elevates the efficacy of anti-PD-1 anti-tumor treatment and the cytotoxic effects of effector T cells due to the metabolic reprogramming [7]. Therefore, targeting appropriate metabolic pathways can stimulate inflammatory immunity and block tumor metabolism, which are crucial for improving immunotherapy, and GLU metabolism is among those potential targets.
In current years, bioinformatics methods have been extensively used for the analysis of microarray and high-throughput data, including the identification of differentially expressed genes (DEGs). These methods are also useful for analyzing the underlying pathways of various human diseases. Since metabolic abnormalities and immune cell infiltration are important in the etiology of PE, the present study used the comprehensive genomic assessment of public databases to find new early diagnostic biomarkers and therapeutic targets by identifying key GLU-DEG modules (genes related to GLU metabolism), pathways, and immune infiltrates to provide ideas for promoting further research in this field.
Materials and methods
Data sources and preprocessing
This investigation used publicly available data, mainly imported from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih. gov/geo/). The flowchart of this study can be found in S1 Fig. Genome-wide expression data related to PE were downloaded using the R package ‘GEOquery’ from the GEO database. The GSE186257 dataset comprised 26 PE patients and 18 controls, while the GSE96984 dataset included 3 PE patients and 4 controls, GSE24129 contained 8 PE patients and 8 controls, and GSE35574 included 19 PE patients and 40 controls. Non-biological technical biases inducing batch effects were corrected with the help of the ComBat protocol of the “sva” R package [8]. Moreover, the correction degree was assessed using principal component analysis (PCA). This investigation honored each database’s data access policies.
A total of 134 GLU metabolism-linked genes were acquired from the Msigdb database (http://www.gseamsigdb.org/gsea/msigdb/index.jsp) [9,10].
PE-associated DEGs
The DEGs between control (n = 22) and PE (n = 29) samples were identified using the “limma" package (version 3.50.0) in R [11]. The parameters were set as |log2 fold change| > 0.5 and adjusted p < 0.05. A heatmap was established using the “pheatmap” package in R with a complete linkage clustering method and Euclidean distance.
Gene Set Variation Analysis (GSVA)
The GSVA is a non-parametric, unsupervised gene enrichment method that utilizes gene expression profiles to identify the correlation of genetic features with biological pathways. To elucidate the difference in the biological role between the control and PE cohorts, GSVA was carried out with “c2.cp.kegg.v7.5.1.symbols” using the R “GSVA" (v 1.42.0) package. The “pheatmap" (v 1.0.12) package in R was used for visualization of the results. Furthermore, 50 HALLMARK gene sets were imported from the MSigDB database (http://software.broadinstitute.org/gsea/msigdb) and used as a reference. The ssGSEA function in GSVA was utilized to assess the GSVA score of each gene set in various samples. These scores were then compared between the control and PE groups using the "limma" package.
Weighted Gene Co‑expression Network Analysis (WGCNA) and assessment of significant modules
The WGCNA algorithm in the R WGCNA package (v 1.70–3) was employed for the construction of co-expression networks [12]. Similarities in gene expressed were analyzed using Pearson correlation coefficients and were weighted using a power function to produce a scale-free network. The ‘PickSoftThreshold’ package in R was used to establish a weighted adjacency matrix by increasing the co-expression similarity to a power β = 22. In terms of co-expression, a gene module represents a densely interconnected gene cluster. WGCNA utilizes hierarchical clustering to assess gene modules, which are indicated by different colors. Furthermore, various modules were identified with the help of the dynamic tree-cut method. Cluster analysis was utilized for module selection, during which the adjacency matrix (a topology similarity measure) was transformed into a topology overlay matrix (TOM). To identify the link between modules and GLU, the relationship of the module eigengene (ME, the module’s 1st principal element that depicts its overall expression level) to GLU was elucidated by Pearson’s correlations. Modules showing significant links with GLU were identified, and heatmaps of the topological overlap in the gene network were generated to observe co-expression module structure. The hierarchical clustering dendrogram of the eigengenes and the corresponding eigengene network heatmap were used to summarize the links between modules. The GLU-related DEGs were acquired by intersecting DEGs and genes from the GLU-related module.
Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment and Gene Ontology (GO) analyses
The GO categories are molecular function (MF), biological process (BP), and cellular component (CC). KEGG is a bioinformatics tool for identifying pathways in which genes are significantly enriched [13]. The KEGG enrichment and GO analyses were performed using the “clusterProfiler" (v 4.2.2) package in R (p-value < 0.05) using the GLU-related DEGs [14].
GeneMANIA
The GeneMANIA website (http://genemania.org) was used to construct protein-protein interaction (PPI) networks. GeneMANIA can also predict the link between hub genes and functionally similar genes, such as protein-DNA interaction or PPI, biochemical and physiological reactions, pathways, co-expression, and co-localization [15].
Receiver Operating Characteristic (ROC) curve
The ROC curve is a plot with test sensitivity on the y-coordinate vs. its 1-specificity or false positive rate (FPR) on the x-coordinate. ROC curves can effectively elucidate the efficiency of diagnostic assays. The most common metric acquired from the ROC plot of sensitivity vs. 1–1-specificity is the area under the curve (AUC). Here, “pROC” [PMID: 21414208] in R was utilized to generate ROC curves to elucidate the AUC values for identifying signature genes and elucidating their diagnostic potential [16]. It was assessed on a 0.5 (a “coin flip”) to 1 (perfect discrimination) scale. Generally, an AUC of > 0.9 = outstanding, 0.8–0.9 = excellent, 0.6–0.8 = acceptable, and 0.5 = no discrimination.
Immune infiltration analysis
The ssGSEA is a variation of the GSEA algorithm that, instead of assessing enrichment scores for various groups (disease vs. control) and gene sets (i.e., pathways), provides a score for each gene set pair and sample [17]. Each ssGSEA enrichment score indicates how many genes in a specific set within a sample are down or up-regulated.
Based on the 28 types of immune cells acquired from the Tumor and Immune System Interactions Database (TISIDB) (http://cis.hku.hk/TISIDB/index.php) [18], including eosinophils, regulatory T cells, activated, immature, and memory B cells, natural killer (NK) cells, gamma delta T cells, CD56 bright and dim NK cells, Type 1, 17, and 2 T helper cells, myeloid-derived suppressor cells, central memory, activated, and effector memory CD8 T and CD4 T cells, T follicular helper cells, macrophages, NK T cells, activated, plasmacytoids, and immature dendritic cells, monocytes, neutrophils, and mast cells, the relative enrichment score of each type of immune cell was determined from each sample’s gene expression profile. Variations in the immune cell infiltration levels between PE and control samples were evaluated using the R package "ggplot2" (v 3.3.6) [19].
Construction of the RBP–mRNA network
StarBase (https://starbase.sysu.edu.cn/tutorialAPI.php#RBPTarget) is a widely utilized open-source platform for assessing the interactions of ncRNA via degradome-seq, RNA-RNA interactome, and CLIP-seq data. StarBase was utilized to elucidate the link between the expression of mRNA and RNA-binding proteins (RBPs). A clusterNum ≥ 5, p-value < 0.05, and clipExpNum ≥ 5 were used as the cutoff parameters for determining the key mRNA-RBP pairs in PE. Furthermore, using Cytoscape, the RBP-mRNA network was generated.
Sample collection
Samples were acquired from women who were registered in the Department of Obstetrics and Gynecology of the Gansu Provincial Hospital for cesarean section from April 1 to September 30, 2023. This research was authorized by the Medical Ethics Committee of the Gansu Provincial Hospital (approval no. 2023–039). Participants were first informed about the research, and then provided written consent. There were two cohorts, namely, health control (HC) and PE. The placental tissue was sampled immediately after the surgery, with each sample divided into two portions: one immersed in RNAsolid tissue RNA stabilization reagent (Servicebio,China) and stored at -80°C until utilized for RT-PCR; the other portion was fixed in polyformalin tissue fixative (Servicebio,China), embedded in paraffin, and utilized for immunohistochemistry.
Immunohistochemistry (IHC)
The paraffin sections were de-waxed and rehydrated,and antigen retrieval was carried out by heating in a microwave oven in citrate buffer. Endogenous peroxidase activity was quenched with 3% hydrogen peroxide, and non-specific binding was blocked with 5% goat serum (Solarbio,China). The slides were then incubated overnight at 4°C with primary antibodies targeting ACSS1 (Proteintech,China), SLC27A2 (Thermofisher,USA), ADAMTS19 (Abcam,USA) and LGR5 (Theroma,USA). Following PBS washing steps, secondary antibodies were incubated, and subsequent incubation with freshly prepared DAB (Solarbio,China) working solution ensued. Counterstaining was performed using hematoxylin, and examined under a microscope.
RT-qPCR
Total RNA was extracted from the placental tissue after cesarean section using TRIzol (TaKaRa, Japan). The cDNA was synthesized per the kit’s protocol (TaKaRa). GAPDH was used as the internal reference, and relative expression of genes was calculated using the 2−ΔΔCT method. qRT-PCR was repeated three times for each sample. The primer sequences used are shown below:
GAPDH
Forward: 5′- TCAAGAAGGTGGTGAAGCAGG-3′,
Reverse: 5′- AAAGGTGGAGGAGTGGGTGTC-3′.
ASSC1
Forward: 5′- TGGCTCACAGGACAGACAACAAG-3′,
Reverse: 5′- GGCGGCATAGAGCAGGTAGC-3’.
LGR5
Forward: 5′- CTCGGTGTGCTCCTGTCCTTG-3′,
Reverse: 5′- GAGGTGAAGACGCTGAGGTTGG-3’.
SLC27A2
Forward: 5′-GCTCCTGGTGAACCTCTGCTG-3′,
Reverse: 5′- AAGGCTTGTGTGGCGTCTGG-3’.
ADAMTS19
Forward: 5′- GAATGGACCCCTTGTTCACGA-3′,
Reverse: 5′-TACATCTAACGGTCCGATGACG-3’.
Statistical analysis
Statistical assessments were conducted using R software v4.1.2. Associations between parameters were assessed by Spearman correlations. Differences between two groups were analyzed using Wilcoxon tests, while differences between multiple groups were assessed by Kruskal–Wallis tests. P-values < 0.05 were considered statistically significant.
Results
Weighted gene co-expression network construction and module identification
The identification of GLU-related gene sets was undertaken using WGCNA. The scale independence and mean connectivity analysis indicated that the average degree of connectivity was close to 0, and scale independence was > 0.85 at the weighted value = 22 (Fig 1A). Altogether, 4 co-expressed modules were determined, and a gray module was assigned for uncorrelated genes, which was subsequently ignored (Fig 1B). To investigate the link between the modules and assess their correlation, the MEs were correlated. Fig 1C illustrates the eigengene network via a dendrogram and a heatmap plot. The topological overlap in the gene network was represented in a heatmap plot (Fig 1D). To understand the physiological importance of the modules, 4 MEs were correlated with GLU, and the most significant associations were assessed. The module-trait correlation heatmap (Fig 1E) indicates that the brown module gene cluster (n = 153) had the strongest positive association with GLU (r = 0.41, p < 0.05). Therefore, the brown module was used in subsequent analyses as it indicated GLU more accurately. Fig 1F indicates the scatterplots of gene significance (GS) for the trait of GLU vs. module membership (MM) in the brown module. MM and GS for GLU had markedly positive associations (cor = 0.37, p < 0.05), suggesting that the brown module’s central (most important) elements were highly correlated with the GLU trait.
(A) Scale-free topological fit index (R2) and soft threshold β = 22. (B) Gene expression network analysis in PE determined specific co-expression data modules. (C) Associations between modules. Top: hierarchical clustering of module eigengenes summarizing the modules identified via the clustering test. Branches of the dendrogram (the meta-modules) group positively correlated eigengenes. Bottom: Heatmap plot of the adjacencies in the eigengene network. In the heatmap, each column and row indicates one module eigengene (color labeled). In the heatmap, blue = low adjacency and red = high adjacency. Red colored squares observed along the diagonal represent the meta-modules. (D) Heatmap of topological overlap in the gene network. In the heatmap, each column and row indicates a gene; progressively darker red color indicates higher topological overlap while light color indicates low topological overlap. The darker squares along the diagonal represent modules. The module assignment and gene dendrogram are indicated along the top and left. (E) Association of GLU with consensus module eigengenes. In the table, each column prepresents a sample or trait, and each row represents a consensus module. Numbers in the table indicate the link between the traits and the corresponding module eigengenes. The p-values are provided in parentheses below the correlations. Based on the color legend, the table is color-coded by correlation. (F) Correlation between gene significance (GS) and module membership (MM) for GLU of genes of the brown module. ‘Cor’ depicts the absolute correlation coefficient between MM and GS.
DEG identification
The comparison of the PE samples and controls identified 356 significantt DEGs between the two cohorts (adjusted p-value < 0.05, |Log2-fold change| > 0.5). In PE specimens, 128 and 228 genes were downregulated and upregulated, respectively. A volcano plot was employed for visualization of the DEGs (Fig 2A). The top five upregulated (BCL6, SPAG4, TMEM45A, QPCT, and HTRA4) and downregulated DEGs (ERICH5, NAV2-AS4, KCTD17, DUSP4, and RAB15) are shown in the heatmap (Fig 2B). Wilcoxon tests showed that the expression of these top 10 genes differed significantly between the two cohorts (p < 0.05, Fig 2C).
(A) Volcano plot of DEGs between the control and PE samples. Gray dots indicate downregulated, upregulated, and insignificant gene expression levels. (B) Heatmap of the top 5 downregulated and upregulated DEGs, respectively. (C) Differences in the expression of the top 10 DEGs between the control and PE cohorts, determined by Wilcoxon tests. Asterisks indicate p-values (*p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001).
A total of 7 GLU-related DEGs were identifed in the intersection between DEGs and GLU-related module genes; these were considered hub genes.
GSVA
To further explore the functional annotations between control and PE samples, GSVA analyses were performed to determine differences in pathways between the two cohorts. Multiple differentially expressed pathways were identified and then visualized using heatmaps. Compared with the control groups, the pathways associated with KEGG_BETA_ALANINE_METABOLISM and KEGG_CARDIAC_MUSCLE_CONTRACTION were markedly reduced in the PE group, whereas the expression of the KEGG_RIG_I_LIKE_RECEPTOR_SIGNALING_PATHWAY and KEGG_ADIPOCYTOKINE_SIGNALING_PATHWAY pathways was notably higher (Fig 3).
Enrichment analyses (GO/KEGG)
To determine the biological activities of the GLU-related DEGs, GO and KEGG enrichment analyses were performed. The GO analysis revealed that the genes were strongly enriched in the BP category sulfur compound transport (GO:0072348) as well asacid-thiol ligase activity (GO:0016878), formation of carbon-sulfur bonds (GO:0016877), ligase activity, and CoA-ligase activity (GO:0016405) in the MF category (Fig 4A). Furthermore, the enriched KEGG pathways included glyoxylate and dicarboxylate metabolism (hsa00630), propanoate metabolism (hsa00640), and pyruvate metabolism (hsa00620) (Fig 4B).
(A) GO enrichment. (B) KEGG enrichment.
Validation of the hub genes
The ROC analysis was utilized to further confirm the diagnostic value of the hub genes. The results showed that ILDR1 (AUC = 0.759), SLC13A3 (AUC = 0.754), ADAMTS19 (AUC = 0.74), SLC27A2 (AUC = 0.73), ACSS1 (AUC = 0.724), LGR5 (AUC = 0.705), SLC13A4 (AUC = 0.702) had similar AUC values (Fig 5A–5G), suggesting that the selected hub genes had the ability to discriminate between the samples and could thus be used as potential biomarkers for PE.
(A) ILDR1, (B) SLC13A3, (C) ADAMTS19, (D) SLC27A2, (E) ACSS1, (F) LGR5, and (G) SLC13A4.
Assessment of interactions between proteins
The GeneMANIA database was used to generate a PPI network of the proteins encoded by the hub genes, identifying seven genes (Fig 6A). Furthermore, to elucidate the role of these genes, KEGG and GO were used to analyze enrichment of 27 genes, including 7 hub genes and 20 related genes. The KEGG analysis showed that the most significantly enriched pathways included those associated with butanoate metabolism (hsa00650), glyoxylate and dicarboxylate metabolism (hsa00630), the PPAR signaling pathway (hsa03320), insulin resistance (hsa04931), propanoate (hsa00640), fatty acid biosynthesis (hsa00061), pyruvate metabolism (hsa00620), valine, leucine, and isoleucine degradation (hsa00280), and fatty acid metabolism (hsa01212) (Fig 6B). The GO data indicated the strong enrichment of these genes in the biological processes of phosphatidylethanolamine biosynthetic process (GO:0006646), positive regulation of hormone secretion (GO:0046887), and heat generation (GO:0031652). In the CC category, enrichment was observed in peroxisomal membrane (GO:0005778), basal plasma membrane (GO:0009925), phosphatidylglycerol biosynthetic process (GO:0006655), microbody membrane (GO:0031903), and basal part of the cell (GO:0045178) and in the MF category, vitamin transmembrane transporter activity (GO:0090482), protein-hormone receptor activity (GO:0016500), and modified amino acid transmembrane transporter activity (GO:0072349) (MF) (Fig 6C).
(A) Gene co-expression network characterization. (B) KEGG and (C) GO analyses of co-expressed genes.
Validation of the gene set
The enrichment of the gene sets carbon metabolism, sulfur compound transport, and PPAR signaling pathway was validated by measurement of the expression levels of the key genes. This showed that the expression of most key genes was reduced in the PE cohort (Fig 7A–7C).
(A) Carbon metabolism, (B) sulfur compound transport, and (C) PPAR signaling in PE and control samples.
Immune cell infiltration
The infiltration of immune cells plays a critical role in PE pathogenesis. Therefore, the associations of control and PE samples with immune cell infiltration were investigated. The PE group indicated a substantially increased infiltration of two types of immune cells, namely, activated B cells and T follicular helper (TFH) cells compared with the control group (Fig 8A). As shown in Fig 8B, the overall infiltration levels of immune cells varied substantially between the PE and control groups. The infiltration of five of the overall 28 immune cell types differed significantly between the two cohorts (p < 0.05) (Fig 8B).
(A) Estimated proportions of infiltrating immune cells in the PE and control cohorts. (B) Heatmap indicating changes in immune cell infiltration between the PE and control cohorts. (C) Correlations between the ACSS1 and TFH cells. (D) Correlations between ACSS1 and memory B cells. (E) Correlations between immune cells. Asterisks represented p-values: *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
Furthermore, associations between the individual hub genes and immune cell infiltration were also assessed (Fig 8C and 8D). It should be noted that, ACSS1 was significantly associated with TFH cells (R = -0.512, p < 0.001) (Fig 8C), and memory B cells (R = 0.52, p < 0.001) (Fig 8D).
Subsequently, correlations between the various type of infiltrating immune cells were estimated. The majority of the immune cells were positively correlated with one another (Fig 8E).
Signaling pathways associated with signature genes
With the help of GSVA, the variability between PE and controls was assessed in 50 HALLMARK signaling pathways. In PE patients, one HALLMARK signaling pathway was found to be significantly up-regulated, namely, HALLMARK_HYPOXIA, and no pathways were down-regulated (Fig 9A).
(A) Comparison of the 50 HALLMARK signaling pathways between the PE and control cohorts. (B) Correlation of the 50 HALLMARK signaling pathways with the target genes. *p < 0.05, **p < 0.01, ***p < 0.001, ****p < 0.0001.
Furthermore, the correlations between 50 HALLMARK signaling pathways and the top five differentially expressed hub genes were elucidated. This showed that ACSS1 was associated with various pathways, including HALLMARK_E2F_TARGETS (Fig 9B).
Construction and functional annotation of crosstalk between RBPs and hub mRNAs
As RBPs interact with mRNAs, 7 hub mRNAs were selected using StarBase, and 7 mRNA/RBP pairs were identified and imported. An RBP-mRNA network was established on the basis of the link between target genes acquired from the online dataset. This network comprised 7 mRNAs, 40 RBPs, 47 nodes, and 105 edges. Fig 10 illustrates the network.
Pink = mRNAs and blue = RBPs.
Verification of hub genes
The value of the hub genes as potential biomarkers was assessed using ROC curve analysis of two independent GEO datasets, GSE24129 and GSE35574, to validate their diagnostic performance in PE. The AUC values that in the dataset GSE35574 (Fig 11A–11D), LGR5 (AUC = 0.75), SLC27A2 (AUC = 0.691), ACSS1 (AUC = 0.679), ADAMTS19 (AUC = 0.65), and in the dataset GSE24129 (Fig 11E–11H), ACSS1 (AUC = 0.844), LGR5 (AUC = 0.828), SLC27A2 (AUC = 0.719), ADAMTS19 (AUC = 0.688), indicated that that the hub genes ADAMTS19, SLC27A2, ACSS1, and LGR5 were potential biomarkers of PE.
LGR5. (A), SLC27A2 (B), ACSS1 (C), ADAMTS19 (D) in the GSE35574 dataset and ACSS1 (E), LGR5 (F), SLC27A2 (G), ADAMTS19 (H) in the GSE24129 dataset.
Verification of hub genes in clinical samples
The expression levels of the hub genes in the placental tissue of patients in the PE and HC groups were determined by quantitative RT-PCR and Immunohistochemical (IHC). Thus indicated that the expression of ADAMTS19, SLC27A2, ACSS 1, and LGR 5 were reduced in PE patients relative to the control group (Fig 12A and 12B), consistent with the bioinformatics analysis. The clinical data of HC pregnant women and PE patients are shown in S1 Table.
(A) Relative mRNA expression of hub genes in the placental tissue of HC and PE. (B) Immunohistochemical analyses of paraffin-embedded placental tissue of HC and PE.****p < 0.0001, **p < 0.01, *p < 0.05.
Discussion
Preeclampsia is a placenta-associated disorder that causes serious complications to both the mother and child, including eclampsia, cerebrovascular accidents, fetal growth restriction, and placental abruption. Due to the absence of early diagnostic markers, women with PE often miss important treatment opportunities and times, resulting in poor prognosis [20,21]. Furthermore, the literature indicates that immune cell infiltration is closely linked with the incidence of PE. A serum metabolomic analysis found clear differences in GLU metabolism between women with normal pregnancies and those with PE [22]. Therefore, identifying specific diagnostic indices and elucidating immune cell infiltration in PE could profoundly improve its prognosis. Bioinformatics is a powerful tool used to identify molecular indices. The present study identified diagnostic indices for PE and explored the role of infiltrating immune cells in PE.
In addition to maintaining the cellular redox state and antioxidant function, glutaminolysis is associated with the biosynthesis of many essential molecules by providing a carbon source for synthetic lipid compounds. Fatty acids are an important substrate for energy supply; therefore, sufficient fatty acid metabolism is essential during pregnancy, and reduced expression of fatty acid metabolic genes in the placenta has been associated with PE [23]. Among the 7 hub genes (ILDR1, SLC13A3, ADAMTS19, SLC27A2, ACSS1, LGR5, and SLC13A4) identified in this study, SLC13A3 is a transporter of citrate, and tricarboxylate, ACSS 1 is a mitochondrial matrix enzyme that promotes acetyl-CoA production, a precursor for fatty acid synthesis [24,25], while SLC27A2 transports fatty acids in oxidative stress and damage to endothelial cells [26]. In addition, citrate is also predictive of PE [27]. Immunoglobulin-like domain-containing receptor 1 (ILDR 1) is a member of the lipoprotein residual receptor family and it has been found that the loss of ILDR 1 results in the inhibition of glucose metabolism in mice [28]. SLC13A4 is a sodium sulfate transporter protein that is mainly expressed in the placenta, and its reduced activity results in defective placental development [29]. ADAMTS19 is a member of the metalloproteinase family that can cleave the components of the extracellular matrix or interact with modulatory factors, affecting cell migration, adhesion, growth, and angiogenesis [30]; an early absence of placental trophoblast invasion is one of the mechanisms associated with PE. Moreover, it has been suggested that ADAMTS19 also affects the biological activity of trophoblast cells; however, further experimental research is required.
The GO and KEGG enrichment analyses identified pathways associated with sulfur compound transport, PPAR signaling, and energy metabolism. Sulfur compounds include sulfate and sulfur-containing amino acids, which are obligate nutrients for healthy development and growth and can be acquired from intracellular sulfur amino acids metabolism and the diet [31]. In PE, the levels of sulfite, sulfur amino acid ergotine, and hydrogen sulfide are markedly reduced [32–34]. Immune functions and metabolism are the core mechanisms used by the body to maintain tissue homeostasis [35]. The metabolism of immune cells is closely related to their phenotype and function [36]. The phenotypic and functional abnormalities of decidual immune cells, such as macrophages and NK cells, can cause dysfunctional trophoblast invasion and defective spiral artery reconstruction, which is the main pathogenesis of PE [37]. PPARs belong to the nuclear hormone receptor family and are implicated in lipid metabolism, immune function, cell growth, and differentiation. Furthermore, free fatty acids utilize PPARs to influence macrophage function and differentiation [38]. Lipids can participate indirectly as substrates or directly as disease mediators in PE etiology [39]. Lipid peroxides can activate endothelial cells and reinforce the pathology of diabetes mellitus and essential hypertension, which are significant risk factors for PE [40]. PPARs can also directly upregulate the expression of fatty acid-binding protein (FABP) and promote fatty acid uptake and accumulation in trophoblast cells, thus inhibiting their invasion ability [41].
GSVA provides valuable data on numerous genes with a relatively smaller fold changes. Here, GSEA analysis was conducted on the gene profiles of the datasets, which revealed multiple gene sets that were markedly enriched in the PE cohort. Among them, the expression of an adipocytokine signaling pathway was significantly increased; this included the synthesis and secretion of many hormones and cytokines by adipocytes to regulate energy metabolism, inflammation, the immune response, cell proliferation, and other physiological functions by binding to cell-surface receptors. The association of various adipokines with insulin resistance and atherosclerosis highlights its importance in metabolic diseases. There are many studies on adipokines and PE. Furthermore, C1q-TNF-related protein 9 (CTRP 9) can affect trophoblast cell proliferation and migration and cause PE [42]. Moreover, Apelin improves endothelial dysfunction to alleviate acute liver and kidney injury in PE [43]. Additionally, a birth credit analysis suggested that adipokine leptin may be a biomarker of PE [44].
To further elucidate the activities of infiltrating immune cells in PE, ssGSEA was employed, which indicated that enhanced infiltration of activated B and TFH cells and reduced infiltration of CD56 bright NK cells, monocytes, and memory B cells might be linked with the incidence of PE. The correlation of 7 hub genes with infiltrating immune cells showed significant negative associations with activated B and TFH cells. Both T and B cells participate in humoral immunity, and increased levels of pro-inflammatory T cells can lead to oxidative stress and cause chronic inflammation characterized by increased production of pro-inflammatory factors and autoantibodies [45]. Women with PE often show superficial chronic inflammatory trophoblast invasion, and the study found that B cell levels were markedly elevated in the peripheral blood of PE patients [46,47]; B cells can induce the production of pro-inflammatory T cells [48]. TFH cells can secrete cytokines to stimulate B cell proliferation. Therefore, it is suggested that the reduced expression of these 7 hub genes might increase the production of activated B and TFH cells and participate in the occurrence and development of PE. These hypotheses require more research to determine the complex communication between genes and immune cells.However, we recognize the limitations in immediate clinical application and the importance of exploring peripheral blood-based biomarkers for practical diagnostic purposes. We will consider expanding our analysis to include peripheral blood samples in future research to enhance clinical relevance.
Conclusion
In summary, this study identified infiltrating immune cells, WGCNA modules, DEGs, enriched pathways, and hub genes that might be associated with PE pathogenesis. These data furnish novel evidence for understanding PE pathogenesis and lays the foundation for future studies aimed at developing clinically relevant diagnostic and therapeutic strategies.
Supporting information
S1 Table. Clinical characteristics of preeclamptic and health control pregnancies.
https://doi.org/10.1371/journal.pone.0303471.s002
(XLSX)
Acknowledgments
We thank all the reviewers who participated in the review and MJEditor (www.mjeditor.com) for its linguistic assistance during the preparation of this manuscript.
References
- 1. Jung E, Romero R, Yeo L, Gomez-Lopez N, Chaemsaithong P, Jaovisidha A, et al. The etiology of preeclampsia. Am J Obstet Gynecol. 2022;226(2S):S844–S866. pmid:35177222
- 2. Tanner MS, Davey MA, Mol BW, Rolnik DL. The evolution of the diagnostic criteria of preeclampsia-eclampsia. Am J Obstet Gynecol. 2022;226(2S):S835–S843. pmid:35177221
- 3. Lu HQ, Hu R. The role of immunity in the pathogenesis and development of pre-eclampsia. Scand J Immunol. 2019;90(5):e12756. pmid:30739345
- 4. Mor G, Cardenas I. The immune system in pregnancy: A unique complexity. Am J Reprod Immunol. 2010;63(6):425–433. pmid:20367629
- 5. Hu XB, Ma ZF, Xu BB, Li SL, Yao ZQ, Liang BS, et al. Glutamine metabolic microenvironment drives M2 macrophage polarization to mediate trastuzumab resistance in HER2-positive gastric cancer. Cancer Commun (Lond). 2023;43(8):909–937. pmid:37434399
- 6. Li W, Ling ZH, Wang JM, Su ZC, Lu JY, Yang X, et al. ASCT2-mediated glutamine uptake promotes Th1 differentiation via ROS-EGR1-PAC1 pathway in oral lichen planus. Biochem Pharmacol. 2023;216:115767. pmid:37634599
- 7. Huang RH, Wang H, Hong J, Wu JY, Huang O, He JR, et al. Targeting glutamine metabolic reprogramming of SLC7A5 enhances the efficacy of anti-PD-1 in triple-negative breast cancer. Front Immunol. 2023;14:1251643. pmid:37731509
- 8. Leek JT, Johnson WE, Parker HS, Jaffe AE, Storey JD. The sva package for removing batch effects and other unwanted variation in high-throughput experiments. Bioinformatics. 2012;28(6):882–883. pmid:22257669
- 9. Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics. 2011;27(12):1739–1740. pmid:21546393
- 10. Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: A knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–15550. pmid:16199517
- 11. Ritchie ME, Phipson B, Wu D, Hu YF, Law CW, Shi W, et al. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015;43(7):e47. pmid:25605792
- 12. Langfelder P, Horvath S. WGCNA: An R package for weighted correlation network analysis. BMC Bioinformatics. 2008;9:559. pmid:19114008
- 13. Kanehisa M, Goto S. KEGG: Kyoto encyclopedia of genes and genomes. Nucleic Acids Res. 2000;28(1):27–30. pmid:10592173
- 14. Yu GC, Wang LG, Han YY, He QY. clusterProfiler: An R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–287. pmid:22455463
- 15. Warde-Farley D, Donaldson SL, Comes O, Zuberi K, Badrawi R, Chao P, et al. The GeneMANIA prediction server: Biological network integration for gene prioritization and predicting gene function. Nucleic Acids Res. 2010;38(Web Server issue):W214–W220. pmid:20576703
- 16. Robin X, Turck N, Hainard A, Tiberti N, Lisacek F, Sanchez JC, et al. pROC: An open-source package for R and S+ to analyze and compare ROC curves. BMC Bioinformatics. 2011;12:77. pmid:21414208
- 17. Wu S, Lv XK, Li YL, Gao XY, Ma ZQ, Fu X, et al. Integrated machine learning and single-sample gene set enrichment analysis identifies a TGF-beta signaling pathway derived score in headneck squamous cell carcinoma. J Oncol. 2022;2022:3140263. pmid:36090900
- 18. Ru BB, Wong CN, Tong Y, Zhong JY, Zhong SSW, Wu WC, et al. TISIDB: An integrated repository portal for tumor-immune system interactions. Bioinformatics. 2019;35(20):4200–4202. pmid:30903160
- 19. Ito K, Murphy D. Application of ggplot2 to pharmacometric graphics. CPT Pharmacometrics Syst Pharmacol. 2013;2(10):e79. pmid:24132163
- 20. Liu YJ, Xu BR, Fan CF. Single-cell RNA sequencing and microarray analysis reveal the role of lipid-metabolism-related genes and cellular immune infiltration in pre-eclampsia and identify novel biomarkers for pre-eclampsia. Biomedicines. 2023;11(8):2328. pmid:37626824
- 21. Chen SR, Ke YM, Chen WH, Wu SJ, Zhuang XX, Lin QY, et al. Association of the LEP gene with immune infiltration as a diagnostic biomarker in preeclampsia. Front Mol Biosci. 2023;10:1209144. pmid:37635936
- 22. Liu XF, Li MD, Lu JJ, Li Y, Zeng AR, Qiang JW. Magnetic resonance spectroscopy and liquid chromatography-mass spectrometry metabolomics study may differentiate pre-eclampsia from gestational hypertension. Eur Radiol. 2023;33(7):4554–4563. pmid:36809432
- 23. Abascal-Saiz A, Fuente-Luelmo E, Haro M, Fioravantti V, Antolín E, Ramos-Álvarez MP, et al. Decreased fatty acid oxidation gene expression in pre-eclampsia according to the onset and presence of intrauterine growth restriction. Nutrients. 2023;15(18):3877. pmid:37764661
- 24. Colas C, Pajor AM, Schlessinger A. Structure-based identification of inhibitors for the SLC13 family of Na+/dicarboxylate cotransporters. Biochemistry. 2015;54(31):4900–4908. pmid:26176240
- 25. Gao X, Lin SH, Ren F, Li JT, Chen JJ, Yao CB, et al. Acetate functions as an epigenetic metabolite to promote lipid synthesis under hypoxia. Nat Commun. 2016;7:11960. pmid:27357947
- 26. Ding XY, Yang Z, Han YW, Yu H. Correlation of long-chain fatty acid oxidation with oxidative stress and inflammation in pre-eclampsia-like mouse models. Placenta. 2015;36(12):1442–1449. pmid:26582505
- 27. Mayrink J, Leite DF, Nobrega GM, Costa ML, Cecatti JG. Prediction of pregnancy-related hypertensive disorders using metabolomics: A systematic review. BMJ Open. 2022;12(4):e054697. pmid:35470187
- 28. Chandra R, Aryal DK, Douros JD, Shahid R, Davis SJ, Campbell JE, et al. Ildr1 gene deletion protects against diet-induced obesity and hyperglycemia. PLoS One. 2022;17(6):e0270329. pmid:35749484
- 29. Rakoczy J, Zhang Z, Bowling FG, Dawson PA, Simmons DG. Loss of the sulfate transporter Slc13a4 in placenta causes severe fetal abnormalities and death in mice. Cell Res. 2015;25(11):1273–1276. pmid:26292834
- 30. Cal S, López-Otín C. ADAMTS proteases and cancer. Matrix Biol. 2015;44–46:77–85. pmid:25636539
- 31. Dawson PA, Elliott A, Bowling FG. Sulphate in pregnancy. Nutrients. 2015;7(3):1594–1606. pmid:25746011
- 32. Zheng XG, Wu WB, Zhou Q, Lian YH, Xiang YQ, Zhao XZ. Targeted bisulfite resequencing of differentially methylated cytosines in pre-eclampsia reveals a skewed dynamic balance in the DNA methylation of enhancers. Clin Sci (Lond). 2023;137(3):265–279. pmid:36645190
- 33. Ho KM. Dose-related relationship between ergothioneine concentrations and risk of preeclampsia. Biosci Rep. 2023;43(8):BSR20231076. pmid:37602425
- 34. Du JX, Wang P, Gou Q, Jin S, Xue HM, Li DS, et al. Hydrogen sulfide ameliorated preeclampsia via suppression of toll-like receptor 4-activated inflammation in the rostral ventrolateral medulla of rats. Biomed Pharmacother. 2022;150:113018. pmid:35483194
- 35. Mathis D, Shoelson SE. Immunometabolism: An emerging frontier. Nat Rev Immunol. 2011;11(2):81. pmid:21469396
- 36. Mor G, Aldo P, Alvero AB. The unique immunological and microbial aspects of pregnancy. Nat Rev Immunol. 2017;17(8):469–482. pmid:28627518
- 37. Li ZH, Wang LL, Liu H, Muyayalo KP, Huang XB, Mor G, et al. Galectin-9 alleviates LPS-induced preeclampsia-like impairment in rats via switching decidual macrophage polarization to M2 subtype. Front Immunol. 2019;9:3142. pmid:30687334
- 38. Gaete PV, Nieves-Barreto LD, Guatibonza-García V, Losada-Barragán M, Vargas-Sánchez K, Mendivil CO. Medium-chain fatty acids modify macrophage expression of metabolic and inflammatory genes in a PPAR β/δ-dependent manner. Sci Rep. 2023;13(1):11573. pmid:37463952
- 39. Kenny LC, Broadhurst DI, Dunn W, Brown M, North RA, McCowan L, et al. Robust early pregnancy prediction of later preeclampsia using metabolomic biomarkers. Hypertension. 2010;56(4):741–749. pmid:20837882
- 40. Wojcik-Baszko D, Charkiewicz K, Laudanski P. Role of dyslipidemia in preeclampsia-A review of lipidomic analysis of blood, placenta, syncytiotrophoblast microvesicles and umbilical cord artery from women with preeclampsia. Prostaglandins Other Lipid Mediat. 2018;139:19–23. pmid:30248406
- 41. Singh S, Dhar R, Karmakar S. Fenofibrate mediated activation of PPARα negatively regulates trophoblast invasion. Placenta. 2022;126:140–149. pmid:35803128
- 42. Li L, Gu ZY, Zhang JJ. CTRP9 overexpression attenuates palmitic acid-induced inflammation, apoptosis and impaired migration in HTR8/SVneo cells through AMPK/SREBP1c signaling. Exp Ther Med. 2022;24(1):459. pmid:35747146
- 43. Hamza RZ, Diab AAA, Zahra MH, Asalah AK, Attia MS, Moursi SM. Ameliorative effect of apelin-13 against renal complications in L-NAME-induced preeclampsia in rats. PeerJ. 2021;9:e11110. pmid:33850656
- 44. Wang YJ, Bai XN, Guo X, Gao XL, Chen YY, Li HR, et al. Bioinformatics analysis combined with clinical sample screening reveals that leptin may be a biomarker of preeclampsia. Front Physiol. 2023;13:1031950. pmid:36685185
- 45. Harmon AC, Cornelius DC, Amaral LM, Faulkner JL, Cunningham MW Jr, Wallace K, et al. The role of inflammation in the pathology of preeclampsia. Clin Sci (Lond). 2016;130(6):409–419. pmid:26846579
- 46. Liao AH, Liu LP, Ding WP, Zhang L. Functional changes of human peripheral B-lymphocytes in pre-eclampsia. Am J Reprod Immunol. 2009;61(5):313–321. pmid:19341382
- 47. Zeng B, Kwak-Kim J, Liu Y, Liao AH. Treg cells are negatively correlated with increased memory B cells in pre-eclampsia while maintaining suppressive function on autologous B-cell proliferation. Am J Reprod Immunol. 2013;70(6):454–463. pmid:24033608
- 48. Muzzio DO, Soldati R, Rolle L, Zygmunt M, Zenclussen AC, Jensen F. B-1a B cells regulate T cell differentiation associated with pregnancy disturbances. Front Immunol. 2014;5:6. pmid:24478775