Identification of Immune-Related Genes as Biomarkers for Uremia

Introduction

Uremia is a symptom of the final stage of chronic renal failure, which is clinically characterized by aberrant water, electrolyte, acid, and base balance and increased blood levels of metabolites, such as creatinine and urea.¹ The uremia stage is often associated with a number of secondary conditions and comorbidities of chronic kidney disease (CKD), including, but not limited to: decreased renal function and accumulation of metabolites; pathological changes in the circulatory system, including anemia, platelet dysfunction, and bleeding; neuromuscular dysfunction and cognitive impairment; vascular endothelial dysfunction and vascular disease progression; endocrine and metabolic disorders, such as insulin resistance, gonadal dysfunction, hyperparathyroidism, and CKD-mineral and bone metabolism disorder; congenital and adaptive immunological diseases characterized by inflammation and immune deficiency.^1–3

Patients in uremia stage are associated with poor prognosis of High morbidity, high prevalence of cardiovascular disease, high mortality rates. The main therapeutic strategies for uremia are kidney transplantation or dialysis (including hemodialysis and peritoneal dialysis). However, the expensive treatment is a major economic burden for society and individuals. Thus, uremia is a major challenge to global health. Recently, the number of studies on uremia has markedly increased. However, the pathological mechanisms of uremia have not been elucidated. The elucidation of the molecular mechanism and the identification of targets for early intervention, which are common concerns of interdisciplinary research and challenging limitations, are the major focus of renal disease research.^1–3

The dynamic regulatory effects of the uremic environment on the immune system are known as uremic inflammation and include immune activation and immunosuppression.⁴ The aberrant activation of the innate immune system, especially monocytes, is one of the features of uremic inflammation.^5,6 Patients with end-stage renal disease treated with renal replacement therapy exhibit premature immunological aging of T cells, which may be the cause of uremia-related immune dysfunction.⁷

This study used bioinformatics methods to examine the correlation between immune-related genes (IRGs) and uremia and screen molecular markers that can aid in the diagnosis, elucidation of pathogenesis, and treatment of uremia.

Materials and Methods

Data Source

The data of an independent uremia cohort (GSE37171) were obtained from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/). The GSE37171 dataset comprised data of blood samples from 75 patients with uremia and 40 healthy subjects. The dataset was divided into training (20 healthy and 63 uremic samples) and validation (20 healthy and 12 uremic samples) sets.³ Additionally, the data of 1793 IRGs were retrieved from the ImmPort database (https://immport.niaid.nih.gov) (Supplementary File 1).

Screening of Differentially Expressed Genes (DEGs)

The DEGs between the uremia and control groups in the training set were screened using the limma package (|log fold-change| > 0.5; adj. p < 0.05).⁸ The screening results were visualized using the volcano plot and heatmap, which were plotted using the R language ggplot2 (version 3.3.6) package and pheatmap package (version 1.0.12), respectively.⁹

Identification of Differentially Expressed IRGs (DEIRGs)

The clinical traits of the uremia and control groups in the training set were used to perform weighted gene co-expression network analysis (WGCNA) using the R package WGCNA to obtain key modules and key module genes.¹⁰ The outlier samples in the training set were filtered out using the Hclust function. Next, the soft threshold of data was determined and the adjacency was calculated. Gene modules were obtained using dynamic tree cutting with a minimum of 30 genes in each gene module, and similar modules were merged. The modules related to uremia were selected as key modules and key module genes for subsequent analysis.

DEIRGs were obtained by intersecting the DEGs, key module genes, and 1793 IRGs using the VennDiagram package.¹¹ Functional enrichment analysis of DEIRGs was performed using the R language clusterProfiler (p < 0.05).¹²

Protein-Protein Interaction (PPI) Network

Protein interactions between DEIRGs associated with uremia were characterized using the STRING database (http://string-db.org) (score > 0.4). The topological properties of the PPI networks were analyzed using cytoHubba in the Cytoscape plugin.¹³ The maximal clique centrality (MCC), degree, and stress algorithms were applied to evaluate and select candidate genes.

Identification and Validation of Biomarkers

In the training set, three machine learning models (generalized linear model (GLM), support vector machine-recursive feature elimination (SVM-RFE), and random forest (RF)) were executed using the R package “DALEX” to plot the residual distribution and obtain the best model. The biomarkers were obtained from the best model. The receiving operating characteristic (ROC) curves of the biomarkers were generated using the survROC R package to evaluate the diagnostic value of the biomarkers.¹⁴ Next, the validation set was used to further assess the diagnostic power of the biomarkers. The expression of biomarkers in the training set and validation set was displayed.

Construction and Verification of a Nomogram

The nomogram of biomarkers in the training set was constructed using the R language rms package. The predictive power of the nomogram model was evaluated using clinical impact curves (CICs) and decision curve analysis (DCA).¹⁵

Evaluation of Infiltrating Immune Cell Landscape

The distribution of immune cells in the training set was examined using the CIBERSORT algorithm.¹⁶ The differential distribution of immune cells in the uremia and control samples was compared using the Wilcox test.¹⁶ Next, Spearman correlation analysis was performed to explore the correlation between biomarkers and differentially abundant immune cells.^17,18

Gene Set Enrichment Analyses (GSEA) and Ingenuity Pathway Analysis (IPA)

Functional enrichment analysis of biomarkers in the training set was performed using the GSEA function of the R package clusterProfiler (p < 0.05; |normalized enrichment score|>1).¹⁹ The gene identifiers and corresponding expression values of the DEGs in the training set were uploaded to the IPA software, which was used to analyze the unique pathways of the biomarkers.²⁰

Validation of Expression of Biomarkers

The expression of biomarkers was verified using quantitative real-time polymerase chain reaction (qRT-PCR). We collected blood samples from 5 uremic patients and 5 healthy subjects for validation at the Guangxi International Zhuang Medical Hospital in February 2023. The 5 uremic patients met the following four conditions: (1) Age 18–90 years old, male or female; (2) Meet the diagnostic criteria for chronic kidney disease, with creatinine clearance <10mL/min; (3) Have been on haemodialysis for more than 3 months and are dialysed 3 times a week for 4 hours each time; (4) No serious cognitive impairment, voluntarily participate in this study, and sign the informed consent form and pass the medical ethical review. This study was approved by the ethical committee of Guangxi International Zhuang Medical Hospital (2023–007-01) and also complies with the Declaration of Helsinki. Informed consent was obtained from all subjects. The blood samples were processed to isolate peripheral blood mononuclear cells (PBMCs) using Ficoll solution. Total RNA was isolated from PBMCs using the TRIzol method. The isolated RNA was reverse-transcribed into complementary DNA (cDNA) using SureScript-First-strand-cDNA-synthesis-kit (Servicebio, China). qRT-PCR analysis was performed with Universal Blue SYBR Green qPCR Master Mix. The PCR conditions were as follows: 95 °C for 60s (initial denaturation), followed by 40 cycles of 95 °C for 20s (denaturation) and 55 °C for 20s (annealing). GAPDH was used as the housekeeping gene. The primer sequences are shown in Table 1. The relative expression levels of target genes were quantified using the 2−ΔΔCT method. The expression levels of target genes were normalized to those of GAPDH.

Table 1 The Primer Sequences of Biomarkers

Statistical Analysis

All statistical analyses were performed using the R programming language. Differences were considered significant at p < 0.05 unless specified otherwise.

Results

Identification of DEIRGs Between Uremia and Control Samples

In the training set, 3464 DEGs (2059 downregulated genes and 1405 upregulated genes) were identified between the uremia and control groups (Figure 1A and B). All samples in the training set exhibited good clustering with no outliers (Figure 2A). The soft threshold power of 17 was selected to calculate adjacencies (Figure 2B). The brown module, which was the most relevant model to uremia, was selected as the key module and comprised 5354 key module genes (Figure 2C and D). The intersection of 3464 DEGs with 1793 IRGs and 5354 key module genes revealed 124 DEIRGs (Figure 2E).

Figure 1 Identification of differentially expressed genes (DEGs). (A) Heatmap of DEGs between the uremia and control groups. Red and blue indicate upregulation and downregulation, respectively. (B) Volcano map of DEGs between the uremia and control groups.

Figure 2 Weighted gene co-expression network analysis (WGCNA). (A) Sample cluster dendrogram of 20 control samples and 63 uremia samples based on their expression profile. (B) Analysis of mean connectivity of each β value from 1 to 20. β = 17 was chosen for subsequent analyses as it has the biggest mean connectivity when the scale-free fit index is up to 0.85. (C) Heatmap of the module-trait correlation. The brown module exhibited the highest correlation with uremia. (D) Scatter plot of correlation between genes and traits in the MEbrown module. (E) Venn diagram representation of the intersection of differentially expressed genes (DEGs), immune-related genes from the ImmPort database, and the brown module. In total, 124 core DEGs were identified.

DEIRGs Were Enriched in Immune Response Regulation

Gene Ontology (GO) analysis revealed that DEIRGs were enriched in T cell proliferation, leukocyte proliferation, leukocyte-mediated immunity, receptor-ligand activity, and regulation of cell development (Figure 3A). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis revealed that DEIRGs were enriched in neuroactive ligand-receptor interaction, T cell receptor (TCR) signaling pathway, neutrophil extracellular trap formation, and PI3K-Akt signaling pathway (Figure 3B).

Figure 3 Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. (A) GO functional enrichment analysis of the differentially expressed immune-related genes (DEIRGs) with the top 10 biological process (BP) and molecular function (MF) terms and the top nine cellular component (CC) terms. The horizontal coordinate shows the ratio of the enriched gene count to the total number of genes input, while the vertical coordinate shows the GO terms. (B) The top 20 significantly enriched KEGG pathways of the DEIRGs.

PPI Network for DEIRGs

The PPI network of DEIRGs comprised 99 sections and 285 edges. PDCD1 interacted with ZAP70, while NGF interacted with PDGFRB (Figure 4A). The top 20 genes identified using the MCC, degree, and stress algorithms were intersected to obtain 9 candidate genes (Figure 4B).

Figure 4 Protein-protein interaction (PPI) network. (A) PPI network of the 99 differentially expressed immune-related genes (DEIRGs). (B) Venn diagram showing the intersection of the candidate key genes obtained using the degree, maximal clique centrality (MCC), and stress algorithms.

Identification of Biomarkers of Uremia

The three models were constructed based on 9 candidate genes. Analysis of the residual distribution of samples revealed that RF was the best model (Figure 5A and B). The model genes with importance greater than 40 were selected as biomarkers. Finally, four biomarkers (PDCD1, NGF, PDGFRB, and ZAP70) were obtained (Figure 5C and D). The area under the ROC curve (AUC) values for the biomarkers in the training set were > 0.9 (Figure 5E). The results of ROC curve analysis of the validation set were consistent with those of the training set (AUC > 0.9) (Figure 5F), indicating that the biomarkers exhibited good diagnostic performance. Compared with those in the control group, the ZAP70 expression levels were downregulated and the PDCD1, NGF, and PDGFRB expression levels were upregulated in the uremia group in both training and validation sets (Figure 5G).

Figure 5 Evaluation of the diagnostic value of biomarkers and validation of hub biomarkers. (A) Cumulative residual distribution of samples. (B) Boxplot of sample residuals. (C) The importance of the variables in the random forest (RF) model, generalized linear model (GLM), and support vector machine-recursive feature elimination (SVM-RFE) model. (D) The importance of RF models. (E and F) Receiving operating characteristic (ROC) curves were used to evaluate the ability of the model to distinguish uremia cases from healthy controls. (G) Validation of hub genes in the GSE37171 datasets. The results of validation set analysis were consistent with those of training set analysis. ***p<0.001.

Nomogram Model of Biomarkers

A nomogram based on four biomarkers (PDCD1, NGF, PDGFRB, and ZAP70) was constructed (Figure 6A). DCA and CIC revealed that the nomogram model based on four biomarkers exhibited better predictive performance than individual biomarkers (Figure 6B and C).

Figure 6 Establishment of the nomogram model in the uremia training cohort. (A) Nomogram model of uremia. (B) Decision curve for the nomogram model. (C) Clinical impact curve for the nomogram model.

Biomarker and Immune Cell Correlation Analysis

The infiltration levels of 22 immune cells were examined in samples in the training set (Figure 7A). The proportions of resting natural killer (NK) cells, M1 macrophages, monocytes, plasma cells, M0 macrophages, and naïve B cells were significantly different between the uremia and control groups (Figure 7B). The expression of PDCD1 was significantly and negatively correlated with the proportion of resting NK cells and naïve B cells and was significantly and positively correlated with the proportion of M0 macrophages. Meanwhile, the expression of NGF was significantly and negatively correlated with the proportion of resting NK cells and naïve B cells and was significantly and positively correlated with the proportion of M1 macrophages. The expression of PDGFRB was significantly and negatively correlated with the proportion of naïve B cells. ZAP70 expression was positively correlated with the proportion of resting NK cells, naïve B cells, and regulatory T cells (Figure 7C).

Figure 7 Analysis of the immune microenvironment. (A) Wind rose graph of immune cell content in control (left) and uremia samples (right). (B) Boxplot of the differential immune cell content between the uremia and control groups. *p<0.05; **p<0.01; ***p<0.001. (C) Correlation between key genes (PDCD1, NGF, PDGFRB, and ZAP70) and immune cells.

Functional Enrichment Analysis of Biomarkers

GSEA results for each biomarker are shown in Figure 8. GO analysis suggested that the biomarkers were enriched in negative regulation of translation, regulation of translation, and ribonucleoprotein complex biogenesis (Figure 8A–D). KEGG analysis revealed that the biomarkers were enriched in neuroactive ligand-receptor interaction, herpes simplex virus 1 infection, nucleocytoplasmic transport, RNA degradation, and spliceosome (Figure 8E–H).

Figure 8 Gene set enrichment analysis (GSEA). (A–D) Gene ontology (GO) enrichment analysis of PDCD1, NGF, PDGFRB, and ZAP70 gene set. (E–H) Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of PDCD1, NGF, PDGFRB, and ZAP70 gene set.

IPA revealed that the biomarkers were involved in 30 pathways. Of these, 19 pathways were inhibited, whereas CD274 cancer immunotherapy pathway was activated (Figure 9A, Supplementary File 2). ZAP70 repressed the NK cell signaling pathway (Figure 9B).

Figure 9 Ingenuity pathway analysis. (A) Classical pathway of differentially expressed genes (DEGs). The abscissa is the path name, while the ordinate is the −log (p-value). The yellow line is the default threshold of p = 0.05. Orange and blue colors indicate the activated and inhibited pathways, respectively. The intensities of Orange and blue indicate the degree of activation and inhibition, respectively. (B) The natural killer (NK) cell signaling pathway was ranked first among the signaling pathways with |z-score| >2. IPA predicted that ZAP70 regulates the NK signaling pathway.

Biomarker Expression Validation

qRT-PCR analysis was performed to verify the expression of biomarkers in the uremia and control groups. The expression of PDGFRB was not significantly different between the uremia and control groups. Meanwhile, the expression levels of other biomarkers in the uremia group were higher than those in the control group (Figure 10A–D).

Figure 10 Validation of the expression levels of biomarkers. (A–D) Expression levels of PDCD1, NGF, PDGFRB, and ZAP70 in the uremia and control groups. ns, not significant;** p<0.01; ***p<0.001.

Discussion

Uremia is not an independent disease but a clinical syndrome associated with the deterioration of kidney function, which is a typical manifestation of kidney failure, CKD, and end-stage renal disease.²¹ Inflammation is a key risk factor for the progression of CKD to uremia.²² Additionally, inflammation is a major etiological factor for other diseases that are independent risk factors for CKD. The aberrant regulation of the innate immune system, especially the downregulation and functional decline of T and B lymphocytes, is reported to be one of the characteristics of uremic inflammation.^6,23–27 In patients undergoing hemodialysis, the contact between peripheral blood monocytes and endotoxins in dialysis membranes or dialysate can lead to the activation of monocytes, exacerbating the inflammatory response. In severe cases, the inflammatory response can adversely affect renal function, aggravate the occurrence of complications, such as malnutrition, anemia, atherosclerosis, and amyloidosis, and increase the risk of cardiovascular diseases.⁴

In this study, DEGs between patients with uremia and healthy subjects were intersected with IRGs to obtain DEIRGs associated with uremia. DEIRGs were subjected to functional enrichment and immune infiltration analyses. Various bioinformatics analysis methods and machine learning methods were used to identify four biomarkers (PDCD1, NGF, PDGFRB, and ZAP70). The predictive and diagnostic values of these four biomarkers and their relevance to the immune cell were evaluated.

GO and KEGG enrichment analysis demonstrated that DEIRGs were mostly enriched in immune-related functions and pathways, such as T cell proliferation. The common immune dysfunctions in patients with uremia include T cell dysfunction and immunosuppression. The two types of adaptive responses are antibody responses elicited by B cells and cell-mediated responses elicited by T cells.²⁸ T cell production, proliferation, and antigen specificity are dysregulated in patients with uremia.¹ Several immune cell types, including dendritic cells (DCs), macrophages, and T lymphocytes, are activated in injured renal tissues, suggesting that immune dysregulation may contribute to the development of CKD.^29–31 Some studies have reported that T cell proliferation may improve the prognosis of patients with uremia. The clearance of plasma toxins through dialysis can increase the number and activity of T cells and may be one of the therapeutic effects of dialysis on immune dysfunctions in patients with uremia.^32,33 Metabolic and endocrine disorders caused by renal dysfunction in patients with uremia, as well as the elevated levels of inflammatory factors, promote an imbalance between protein synthesis and degradation pathways, leading to impaired intracellular signal transduction, nucleic acid synthesis, and translation regulation.^1,29,34,35 Thus, uremia affects multiple biological processes. Further studies are needed to elucidate the specific mechanisms and develop novel therapeutic targets for uremia.

Additionally, this study identified four DEIRGs rarely reported in uremia studies. Of these, PDCD1 exhibited the highest correlation coefficient. Previous studies have demonstrated that PDCD1, which is widely expressed in various cell types, including T cells, B cells, DCs, and macrophages, downregulates T cell responses, resulting in immunosuppression events, including T cell apoptosis, T cell failure, T cell anergy, IL-10 production, and regulatory T cell induction.³⁶ Mizuki Kobayashi reported that PDCD1 polymorphisms were associated with severe immune-related adverse events in patients with metastatic renal cell carcinoma who were treated with nivolumab.³⁷ Increased immune activation can lead to immune-mediated tissue damage and autoimmune diseases. Antineutrophil cytoplasmic antibody (ANCA)-associated vasculitis, a systemic autoimmune disease, promotes the progression of rapidly progressive glomerulonephritis. Samy Hakroush reported that in ANCA-associated renal vasculitis, PDCD1 expression was downregulated in the renal tubule interstitium. Tubulointerstitial PDCD1 deficiency is associated with active ANCA-associated renal vasculitis.³⁸ Specific PDCD1 inhibitors promote ANCA-associated renal vasculitis, suggesting a similar mechanistic link between impaired PDCD1 signaling pathways and disease activation.³⁹ Furthermore, PDCD1 and CD274 deficiencies are associated with lupus-like syndrome phenotypes in mice. The expression levels of EGF, PDGF, transferrin, and folate receptors are upregulated in the kidneys of patients with lupus nephritis and lupus nephritis mouse model.^40,41 Wang et al performed WGCNA and immune infiltration analysis and revealed that PDCD1 is involved in T cell activation, monocyte differentiation, cell-cell adhesion regulation, and lymphocyte differentiation, which was consistent with the enrichment results of Zhou et al. The authors identified PDCD1 as a diagnostic biomarker for acute renal rejection.^42,43 PDCD1 is involved in the pathogenesis of various kidney diseases. Uremia is the outcome of various CKDs. We hypothesize that PDCD1 regulates the pathogenesis of uremia by modulating the activation and apoptosis of T and B cells. Zhang³⁶ and Thibult⁴⁴ reported that PDCD1 regulates T and B cells, which is consistent with the results of bioinformatics analysis in this study.³⁶

Previous studies have reported that NGF is involved in the pathogenesis of renal fibrosis. Renal fibrosis, which is characterized by glomerulosclerosis and tubulointerstitial fibrosis, is a common final manifestation of several chronic renal diseases.⁴⁵ Renal disease associated with progressive fibrosis is characterized by elevated renal and serum levels of NGF. Additionally, the downstream signaling pathway of NGF and TGF-b1 mediates most of the key events involved in renal tubular epithelial-interstitial transformation, which has been supported by Vizza et al.⁴⁶

Renal fibrosis is the outcome of the progression from CKD to end-stage renal disease. PDGFRB is upregulated in mesangial cells, podocytes, vascular smooth muscle cells, renal tubule cells, and interstitial cells in animal kidney disease models and patients with kidney disease.⁴⁷ Focal segmental glomerulosclerosis (FSGS) is a primary podocyte disease that is a secondary consequence of severe or progressive glomerular injury. Specific treatment options are not available for FSGS, and many patients progress to end-stage renal disease. Although podocyte injury was considered to be the cause of FSGS initiation, parietal epithelial cells (PECs) are reported to be one of the main effector cells in FSGS formation. Based on the involvement of platelet-derived growth factor (PDGF) in the fibrosis process, Jia et al reported that short-term or long-term treatment of anti-PDGFRB neutralizing antibodies improved renal function, reduced the proliferation of FSGS and PEC, and promoted fibrosis activation in FSGS mice. These findings confirm that PDGFRB plays an important role in the formation and progression of FSGS.⁴⁸

Most patients with uremia exhibit varying degrees of renal fibrosis, which is one of the main pathological and physiological foundations of uremia.^1–3 Therefore, NGF and PDGFRB may be involved in the occurrence and development of uremia by influencing the process of renal fibrosis. PDGFRB is associated with cell development, proliferation, and transformation. The downregulation of PDGFRB is correlated with increased cell apoptosis. Therefore, premature immune aging of T cells⁷ in patients with uremia may be related to PDGFRB although the specific mechanism has not been elucidated.

ZAP70 is reported to be involved in the pathogenesis of clear cell renal cell carcinoma. Naghdibadi et al used WGCNA to cluster DEGs into eight distinct modules. Clear cell renal cell carcinoma-associated DEGs were mainly enriched in the translocation of ZAP70 to immune synapses, endosomal/vacuolar pathways, vascular wall cell surface interactions, and immune-related pathways.⁴⁹ ZAP70, which is mainly expressed in T cells and NK cells, regulates the TCR signaling pathway. The activation of T cells promotes the CD3 (transmembrane protein)-mediated recruitment of ZAP70 to the phosphorylation domain of the TCR complex, leading to T cell activation and proliferation.⁵⁰ Therefore, we hypothesized that ZAP70 promotes the development of uremia by regulating T cell activation and proliferation.

In this study, clinical samples were collected for experimental verification. Compared with those in healthy participants, the serum levels of PDCD1, NGF, and ZAP70 were significantly upregulated and the serum levels of PDGFRB were downregulated in patients with uremia. This is consistent with the results of bioinformatics analysis, indicating the reliability of the results of this study.

However, this study has some limitations. Firstly, our study was based on a limited number of samples from public databases, and expanding the sample size is a pressing issue. In addition, the effect of sample type on gene expression levels needs to be further ruled out by collecting multiple types of samples. Although we validated the gene expression levels by qRT-PCR, further validation of the analyses through experiments is still necessary. More importantly, according to our analysis results, immune cell-related genes and differential immune cells play an important role in the development and progression of uremia, but the specific mechanism of action needs to be further investigated. This will also be the focus of our subsequent research. In conclusion, we will continue to pay attention to the role of the four biomarkers in the developmental process of uremia and explore their mechanisms of action.

Conclusion

This study identified four IRGs associated with uremia using bioinformatics analysis and machine learning methods. The expression of DEGs between healthy individuals and patients with uremia was verified using clinical samples. These genes may affect the occurrence and prognosis of uremia by regulating immune function. The findings of this study improved our understanding of the mechanisms underlying the occurrence and development of uremia and may aid in developing novel therapeutic strategies for uremia.