Identification of Genes Associated with Decreasing Abundance of Monocytes in Long-Term Peritoneal Dialysis Patients

Introduction

The peritoneum, the largest biological internal membrane, covers the surface of abdominal organs and the walls of the abdominal cavity. Currently, peritoneal dialysis (PD) is considered as an effective method of renal replacement therapy for patients with end-stage renal disease. Compared to hemodialysis, PD offers some advantages such as better preservation of residual renal function and similar survival.¹ However, long-term exposure of the peritoneal membrane to non-biocompatible PD fluid results in structural and functional abnormalities of peritoneum, ultimately leading to PD failure.² Clinically, peritonitis is a common complication of PD, which contributes to reduced technique survival in the short and mid-term. In addition, peritoneal fibrosis remains a serious complication of long-term PD (LPD). Moreover, peritoneal fibrosis leads to loss of peritoneal ultrafiltration capacity, with inflammation as a key event of pathogenesis.³

In the peritoneal tissues of patients who undergo an LPD therapy, peritoneal inflammation and fibrosis are regulated by the cytokine network in the celiac immune system.⁴ It is noted that monocytes play an important role in the systemic inflammatory response to severe injury and infection. It has been demonstrated that the altered expression of monocyte subsets affects the modulation of inflammation and immunity with significant clinical implications.⁵ During inflammation, monocytes circulate through the blood and extravasate into inflamed tissues after the general paradigm of the leukocyte recruitment cascade. Upon recruitment to sites of inflammation, monocytes are generally thought to differentiate into macrophages while maintaining the same inflammatory phenotype.⁶ Therefore, it is crucial to fully elucidate the exact pathologic mechanisms underlying the alterations of monocytes and related genes in the peritoneal cells of the peritoneum during LPD, in order to recognize potential therapeutic targets to protect the membrane integrity. In addition, with the development of biotechnology, experimental studies have produced high-throughput transcriptome data and single-cell sequencing data, which play an important role in revealing genes related to diseases or treatments. High-throughput RNA-sequencing (RNA-seq) is an indispensable tool for transcriptome-wide transcriptional variation analysis. RNA-seq can capture various types of RNA, including mRNAs, miRNAs, lncRNAs, and circRNAs.^7,8 High-throughput single-cell transcriptome sequencing, commonly referred to as single-cell RNA sequencing (scRNA-seq), is used to measure gene expression at the single-cell level.⁹ ScRNA-seq analysis of systemic lupus erythematosus nephritis revealed that high expression of interferon-inducible genes in renal tubular cells was associated with disease severity.¹⁰ ScRNA-seq analysis has also been applied to the rejection response after kidney transplantation to deconstruct its heterogeneity and identify new targets for personalized treatment methods.¹¹ In this study, Bulk RNA-seq and ScRNA-seq data from LPD and short-term peritoneal dialysis (SPD) patients were downloaded to identify monocytes-related genes during the PD process.

Materials and Methods

Datasets

In this study, bulk RNA-seq data and scRNA-seq data from peritoneal cell samples from PD patients were used for analysis. GSE125498 dataset (involving RNA-seq data from peritoneal cell samples)¹² and GSE130888 dataset (involving scRNA-seq data from peritoneal cell samples)¹³ were downloaded from the Gene Expression Omnibus (GEO) database. There are 13 long-term PD (LPD) patients and 20 short-term PD (SPD) patients in the GSE125498 dataset. The GSE130888 dataset includes four long-term PD patients, six short-term PD patients, and three normal controls. The GSE125498 dataset was used for differential expression analysis, immune infiltration analysis, and weighted gene co-expression network analysis (WGCNA). The GSE130888 dataset was used to identify monocyte marker genes.

Identification of Differentially Expressed Genes (DEGs) and Immune Infiltration Analysis

The gene expression profile was annotated with the annotation file. DEGs were identified through the limma package in R under the screening criteria of p < 0.05 and |log₂ fold change (FC)| > 0.5. Heat maps and volcano maps of DEGs were, respectively, visualized using “Pheatmap” and “GGPlot” packages in R. In addition, combined with the known immune cell characteristics, the “CIBERSORT” package in R was used to calculate the proportion of 22 kinds of immune cells in the sample based on the expression of the genes. Expression of characteristic genes was extracted to form the expression matrix. The proportion of relevant monocytes was extracted. Correlations between genes and monocytes were calculated by Pearson test and visualized through ggplot packages in R.

Wgcna

WGCNA is a typical systematic biological algorithm to construct gene co-expression networks. The “WGCNA” package in R was used to analyze the co-expression network of genes in the top 25% of the coefficient of variation, followed by the construction of a scale-free gene co-expression network. The “pickSoftThreshold” function was used to select an appropriate soft threshold power regulator to build the scale-free topology. The adjacency matrix is calculated according to the soft threshold. The adjacency matrix is transformed into a topological overlap matrix and the corresponding dissimilarity matrix. The average linkage hierarchical clustering method was used to cluster genes. According to the criteria of the hybrid dynamic shear tree method, the minimum number of genes in each gene network module was set to 30. The characteristic genes of each module were successively calculated. Those hub modules most significantly associated with monocytes were selected. Module connectivity is defined as the absolute value of Pearson correlation between genes. Monocyte relationship is defined as the absolute value of Pearson correlation between each gene and monocyte. In the hub module, genes with module connectivity > 0.8 and monocyte relationship > 0.2 were considered as monocytes related genes.

Identification of Monocytes Marker Genes

ScRNA-seq data in the GSE130888 dataset was analyzed using the “Seurat” and “SingleR” packages in R. A total of 77,197 cells from 10 PD samples were included in the analysis. Quality control was carried out according to the following standards: (1) genes detected in less than 3 cells were excluded; (2) cells with a total number of genes less than 50 were excluded; (3) cells with mitochondrial expression genes more than 5% were excluded; (4) all cells with gene expression less than 1000 were excluded. After this procedure, a total of 66,350 cells were retained. Gene expression in these cells was normalized using linear regression models. The t-distributed random neighborhood embedding algorithm is applied to reduce the dimension of top 20 principal components. Cluster analysis was performed for all cells. Different cell clusters were identified and annotated using a singleR package in R based on the compositional pattern of the marker genes. Differential expression analysis of all genes within cell clusters was performed using the limma package in R to identify marker genes for each cluster. Monocytes marker genes were identified under the threshold value of false discovery rate (FDR) < 0.05 and | log₂FC | > 1.

Identification and Functional Analysis of Core Genes

The DEGs, monocytes related genes in the hub model and monocytes marker genes were intersected to obtain core genes related to monocytes changes during PD. Database for Annotation, Visualization, and Integrated Discovery (DAVID) database was used for functional analysis of core genes. Significantly enriched Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) terms were identified under the screening criteria of p < 0.05 and visualized by bubble maps.

Expression Validation of Core Genes by Reverse Transcriptase Polymerase Chain Reaction (RT-PCR)

RT-PCR was used to validate the expression of core genes in blood samples from LPD and SPD patients. Based on the standards in the GSE125498 and GSE130888 datasets, the inclusion and exclusion criteria for PD patients in this study were set. For LPD patients, the duration of dialysis was more than 25 months. For SPD patients, the duration of dialysis was 0–24 months. The common inclusion criteria of LPD and SPD patients were as follows: (1) patients had no signs of infection at the time of examination and were free of peritonitis for at least 4 weeks prior to the inclusion; (2) patients’ clinical laboratory data were complete. The exclusion criteria of LPD and SPD patients were as follows: (1) patients with peritonitis; (2) patients had severe infectious diseases, systemic inflammatory diseases, autoimmune diseases, and malignant diseases; (3) patients with diabetes; (4) patients had a family history; (5) patients had incomplete clinical data. Based on the above inclusion and exclusion criteria, a total of 14 LPD patients and 11 SPD patients were enrolled from the population examined and treated in General Hospital of Northern Theater Command. The blood samples were collected from July 4, 2022 to November 19, 2022. Clinical information about these patients is listed in Table 1. Blood samples from these patients were collected for RT-PCR. GAPDH and ACTB were used as internal references. 2^−ΔΔct method was used to calculate the expression of gene. 2^−ΔΔct > 1 and 2^−ΔΔct < 1 represent up-regulation and down-regulation, respectively. The study was approved by the ethics committee of the General Hospital of Northern Theater Command. In addition, all individuals provided the informed consent of the patients and their families.

Table 1 Clinical Information of Patients in RT-PCR

Statistical Analysis

T-test was used for statistical analysis in the expression validation in RT-PCR. P < 0.05 was considered as a statistical difference.

Results

Identification of DEGs and Immune Infiltration in LPD Patients

Compared with SPD patients, 823 DEGs were identified in peritoneal cells of LPD patients, including 513 up-regulated and 310 down-regulated genes. Volcanic maps and heat maps of top 100 genes are presented in Figure 1A and B, respectively. In addition, estimation of immune cells in peritoneal cells from PD patients was performed using CIBERSORT package in R. Based on the clustering of 22 kinds of immune cell subsets, peritoneal cell samples were divided into two discrete groups (Figure 2A). It is noted that monocytes were the most abundant immune cell. In addition, the correlation between 22 kinds of immune cell types was further estimated (Figure 2B), which showed strong correlations among some immune cell subpopulations. For example, monocytes were negatively correlated with many T cells. Among which, monocytes were most negatively correlated with T cells CD4 memory resting (r = −0.48, p < 0.05). Moreover, the number of monocytes was significantly decreased in LPD compared with SPD (p = 0.05) (Figure 2C). It is indicated that the number of monocytes decreases with the prolongation of PD time.

Figure 1 Volcanic maps (A) and heat maps (B) of top 100 genes in LPD patients.

Figure 2 Immune infiltration analysis. (A) the heat maps of immune cell subsets; (B) correlation analysis between immune cell types; (C) percentages of immune cell subsets.

Identification of Monocytes Related Genes in the WGCNA

In the top 25% of the coefficient of variation in the sample, 4391 genes were utilized to construct a weighted gene co-expression network. The clustering result showed that there was an outlier sample (GSM3575608) (Figure 3A). After removing outlier sample, the dendrogram of 32 samples and heat map of monocyte abundance are shown in Figure 3B. The β = 6 was chosen as the soft threshold to construct the scale-free network (Figure 3C). Dynamic tree cutting method was utilized to merge modules with dissimilarity < 30%. Totally, 16 modules were finally determined (Figure 3D). It is worth mentioning that genes in the yellow module (involving 263 genes) were the highest negatively correlated with monocytes (R² = −0.8, p < 0.001) (Figure 4A). The scatter plot of genes in yellow module is shown in Figure 4B. In the yellow module, 64 genes were associated with changes in monocytes during PD.

Figure 3 WGCNA in PD patients. (A) sample clustering. Above the red line was an outlier sample; (B) sample dendrogram and heat map of monocyte abundance; (C) filtering of soft thresholds; (D) identification of modules.

Figure 4 Identification of key modules associated with monocytes. (A) correlation between monocytes and different gene modules; (B) scatter plot of genes in yellow module.

Identification of Monocytes Marker Genes

After excluding low-quality cells and genes, a total of 66,350 cells (involving 33,694 genes) from 10 PD patients were used for analysis (Figure 5A). There were 2000 highly variable genes in the variance analysis (Figure 5B). Principal components analysis (PCA) was performed to identify available dimensions and screen for relevant genes. There was no significant segregation between cells (Figure 5C). Top 20 principal components were selected for subsequent analysis (Figure 5D). Through a t-distributed random neighborhood embedding algorithm, the peritoneal cells were successfully divided into 17 independent clusters (Figure 5E). Clusters of 2, 3, 4, 6, 7, 8, 12, 13, 14, and 16 were annotated as monocytes. Interestingly, the number of monocytes in LPD patients was significantly lower than that in SPD patients. The results were consistent with those of the immune infiltration analysis results in the GSE125498 dataset. There were different gene expression profiles in different cell clusters. Some genes were differentially expressed among the 17 clusters. These clusters were annotated by singleR based on the expression patterns of marker genes (Figure 5F). After differential expression analysis between monocytes clusters and other cell clusters, a total of 231 monocytes marker genes were identified.

Figure 5 Identification of monocytes marker genes. (A) 66,350 cells were used for analysis; (B) highly variable genes in the variance analysis; (C) PCA; (D) analysis of top 20 principal components; (E) 17 independent clusters of peritoneal cells; (F) annotation of clusters.

Identification and Functional Analysis of Core Genes in LPD Patients

A total of 16 common up-regulated genes were identified between 823 DEGs, 64 monocytes related genes in the yellow model, and 231 monocytes marker genes (Figure 6A). These genes were granzyme A (GZMA), granzyme M (GZMM), TBC1 domain family member 10C (TBC1D10C), SPARC (osteonectin), cwcv and kazal-like domains proteoglycan 2 (SPOCK2), CD2 molecule (CD2), CD6 molecule (CD6), CD3 epsilon subunit of T-cell receptor complex (CD3E), CD3 delta subunit of T-cell receptor complex (CD3D), CD3 gamma subunit of T-cell receptor complex (CD3G), LCK proto-oncogene, Src family tyrosine kinase (LCK), interferon stimulated exonuclease gene 20 (ISG20), linker for activation of T cells (LAT), ETS proto-oncogene 1, transcription factor (ETS1), spectrin repeat containing nuclear envelope protein 2 (SYNE2), septin 1 (SEPTIN1), and cystatin F (CST7). PPI results showed that there were interactions between these common core genes, except for SYNE2 (Figure 6B). Based on GO analysis of 16 common core genes, T cell activation, immunological synapse, and T cell receptor binding was the most significantly enriched biological process, cytological component, and molecular function, respectively (Figure 6C). In addition, Th1 and Th2 cell differentiation, T cell receptor signaling pathway, Th17 cell differentiation, and hematopoietic cell lineage were significantly enriched signaling pathways (Figure 6D). It is noted that LCK, CD3G, CD3E, CD3D, and LAT were involved in the signaling pathways of Th1 and Th2 cell differentiation, T cell receptor signaling pathway, and Th17 cell differentiation. CD2 was involved in hematopoietic cell lineage signaling pathway.

Figure 6 Identification and functional analysis of core genes. (A) Venn diagram; (B) PPI network; (C) GO analysis; (D) KEGG analysis.

Expression and Correlation Analysis of Core Genes in Monocytes of LPD Patients

In order to explore the expression of 16 core genes in monocytes, relative expression of these genes in each cell cluster was compared. The results showed that the expression of 16 core genes was lower in monocyte clusters than that in other cell clusters (Figure 7). In addition, the correlation between 16 core genes and monocytes was analyzed. The expression of 16 core genes was negatively correlated with the abundance of monocytes (Figure 8).

Figure 7 Expression analysis of core genes in monocytes.

Figure 8 Correlation analysis between core genes and monocytes.

Expression Validation of Core Genes in LPD Patients

Totally, five core genes (CD3D, CD6, CD3E, LAT, and GZMA) were used for expression validation by RT-PCR (Figure 9). The results showed that CD3D, CD6, CD3E, LAT, and GZMA were up-regulated in the blood samples of LPD patients compared with SPD patients. The expression trend consisted with the informatics analysis results.

Figure 9 Expression validation of core genes in LPD patients.

Discussion

In end-stage kidney disease, the altered expression of immune cells including monocyte subsets affects the modulation of immunity and inflammation.¹⁴ Innate immune cells, monocytes, and polymorphonuclear leukocytes recognize pathogens through Toll-like receptors, which trigger phagocytosis, cell activation, and secretion of inflammatory cytokines. A study found that TLR2 and TLR4 expression were significantly up-regulated in monocytes in the end-stage kidney disease group.¹⁵ Furthermore, monocyte count modifies the association between chronic kidney disease and the risk of death.¹⁶ In this study, we found that monocytes were the most abundant immune cell in peritoneal cells of LPD patients. Moreover, the number of monocytes was significantly decreased in LPD compared with SPD. It is indicated that the number of monocytes decreases with the prolongation of PD time. In addition, 16 up-regulated core genes were negatively correlated with the abundance of monocytes, including CST7, SYNE2, SPOCK2, CD2, CD3D, CD6, LCK, ISG20, TBC1D10C, GZMM, CD3E, LAT, GZMA, SEPTIN1, CD3G, and ETS1. Moreover, the expression of 16 core genes was lower in monocyte clusters than that in other cell clusters.

CST7 acts on hematopoietic cell precursors and regulates immune function.^17,18 Significantly increased expression of CST7 has been found in glomeruli.¹⁹ SYNE2 is involved in cell cycle regulation and acted upon by viral proteins.²⁰ SPOCK2 is an important downstream target within the regulation of nephron development.²¹ CD2 is associated with extracellular vesicle characteristics in different kidney diseases.²² In transgenic mice, over expression of RORγt under the control of the CD2 promoter contributes to the development of peritoneum fibrosis.²³ CD3D is up-regulated in the biopsies of acute rejection patients after kidney transplantation.²⁴ CD6 is associated with acute rejection across renal allografts.²⁵ LCK plays an important role in end-stage kidney disease.²⁶ Via acting on T-cell receptor, LCK could be a therapeutic target for acute rejection after kidney transplantation.²⁷ ISG20 is over-expressed in kidney injury and antibody-mediated rejection.^28,29 TBC1D10C is an abundant protein in peripheral blood leukocytes. In immune cells, TBC1D10C plays roles in the regulation of lymphocyte activation.³⁰ Thus, it can be seen that CST7, SYNE2, SPOCK2, CD2, CD3D, CD6, LCK, ISG20, and TBC1D10C play an important role in kidney fibrosis and acute rejection patients after kidney transplantation. However, up to now, there have been few studies on CST7, SYNE2, SPOCK2, CD2, CD3D, CD6, LCK, ISG20, and TBC1D10C regulating the role of monocytes in PD. In this study, we found that the number of monocytes significantly decreased in LPD compared to SPD. The expression of CST7, SYNE2, SPOCK2, CD2, CD3D, CD6, LCK, ISG20, and TBC1D10C were lower in monocyte clusters than that in other cell clusters. Moreover, correlation analysis showed that the expression levels of CST7, SYNE2, SPOCK2, CD2, CD3D, CD6, LCK, ISG20, and TBC1D10C were negatively correlated with monocyte abundance. Therefore, it is hypothesized that CST7, SYNE2, SPOCK2, CD2, CD3D, CD6, LCK, ISG20, and TBC1D10C may play an important role in LPD patients by regulating monocytes, but the specific molecular mechanism needs to be further studied.

GZMM is released in the context of severe inflammation.³¹ GZMM enhances the cytotoxic function and contributes to the immunostimulatory microenvironment.³² CD3E, involved in inflammatory and involved T-cell development, is associated with the transition from acute to chronic kidney injury following ischemia/reperfusion.^33–35 CD3E can predict T-cell–mediated rejection in kidney allograft.³⁶ LAT is involved in inflammatory signaling pathway in chronic PD.³³ GZMA can induce the expression of pro-inflammatory cytokines in macrophages, monocytes, and fibroblasts.^37–41 In kidney, ablation of expression of GZMA leads to an increase in the number of Theileria-containing tumours.⁴² In kidney allografts, GZMA is associated with acute rejection.⁴³ However, there are few studies on GZMM, CD3E, LAT, and GZMA regulating the role of monocytes in PD. In this study, we found that the expression levels of GZMM, CD3E, LAT, and GZMA were up-regulated, and the number of monocytes was significantly decreased in LPD compared to SPD. Moreover, correlation analysis showed that the expression levels of GZMM, CD3E, LAT, and GZMA were significantly negatively correlated with monocyte abundance. Therefore, it is speculated that the increased expression of GZMM, CD3E, LAT, and GZMA may be involved in the peritoneal inflammatory response of LPD, and the mechanism may be related to the decreased abundance of monocytes.

SEPTIN1 is involved in immune and inflammatory processes.⁴⁴ CD3G, associated with inflammation and immune system activation, has been identified as a potential therapeutic target for acute rejection after kidney transplantation.^27,33,45 ETS1, an inflammatory and immune cell-specific transcription factor, plays a key role in kidney development.^46–48 The expression of ETS1 is increased in interstitium and glomeruli during the progression of crescentic glomerulonephritis.⁴⁹ After 12 months of dialysis treatment, CpG sites in ETS1 are significantly less methylated in patients.⁵⁰ In acute kidney injury patients, knocking down ETS1 will alleviate the pyroptosis of renal tubular epithelial cells.⁵¹ In this study, we found that the expression levels of SEPTIN1, CD3G, and ETS1 were up-regulated. Moreover, correlation analysis showed that the expression levels of SEPTIN1, CD3G, and ETS1 were significantly negatively correlated with monocyte abundance. Therefore, it is speculated that SEPTIN1, CD3G, and ETS1 may play a crucial role in both immune and inflammatory processes in LPD by regulating monocyte abundance.

Based on functional analysis of above 16 core genes, we found that LCK, CD3G, CD3E, CD3D, and LAT were involved in the signaling pathways of Th1 and Th2 cell differentiation, T cell receptor signaling pathway and Th17 cell differentiation. CD2 was involved in hematopoietic cell lineage signaling pathway. Zamauskaite et al found that PD patients presented subclinical systemic inflammation, characterized by activated Th1/Th2 responses.⁵² Pharmacological therapeutic strategy that modulates T cell responses could prevent PD fluids-induced peritoneal damage.⁵³ Th17 cells are participated in the pathogenesis of several autoimmune and inflammatory diseases, such as CKD.⁵⁴ Accumulation of Th17 cells in the dialyzed peritoneum is associated with inflammation and fibrosis.⁵⁵ It is found that dialysate cells in PD effluents are largely composed of a small portion of hematopoietic stem cells.⁵⁶ It is indicated that LCK, CD3G, CD3E, CD3D, LAT, and CD2 may be involved in the pathological mechanism of PD process.

In this study, we found that the number of monocytes significantly decreased in LPD compared to SPD, implying that monocytes influence the PD process. Moreover, correlation analysis showed that the expression levels of 16 core genes were negatively correlated with monocyte abundance, which again implied that monocytes play an important role in PD. The negative correlation between these up-regulated core genes and decreased monocytes abundance may be associated with LPD. In addition, functional enrichment analysis showed that the core genes may be involved in the regulation of multiple immune signaling pathways, which laid a theoretical foundation for understanding the molecular mechanism during PD. In short, this study contributes to the understanding of the molecular mechanism of monocytes in LPD and provides potential directions for further research. However, there are limitations to this study. Firstly, the expression of identified genes is needed to be validated in larger numbers of clinical samples. Secondly, a deeper action mechanism between genes, monocytes, and related signaling pathways is needed in further animal models.