Investigating EGF and PAG1 as necroptosis-related biomarkers for diabetic nephropathy: an in silico and in vitro validation study

scRNA-seq data analysis

The Seurat package was utilized for integration and quality control processes, which identified a total of 23,980 cells. Subsequently, the UMAP algorithm revealed twenty-seven distinct clusters, as depicted in Figure 1A. The top 20 most highly variable markers were visualized in Figure 1B. Using gene markers identified from literature relating to kidneys (Supplementary Table 1), we identified 13 cell types from all clusters, as shown in Figure 1C. We then identified 698 down-regulated DEGs and 511 up-regulated DEGs (Supplementary Table 2).

Figure 1. Cell type and the identification of differentially expression genes (DEGs) in the scRNA-seq dataset. (A) The distribution of 27 cell clusters (0-26) in the GSE131882 dataset as determined by the UMAP algorithm. (B) Violin plot showing the expression levels of highly variable gene markers in different cell clusters. (C) Thirteen cell types were identified in the GSE131882 scRNA-seq data. (D) Circle plot of the top 5 GO and KEGG enrichments in the up-regulated DEGs. (E) Circle plot of the top 5 GO and KEGG enrichments in the down-regulated DEGs. Abbreviations: LOH, loop of Henle cell; CD-ICA, collecting duct type A intercalated cell; CT, connecting tubule cell; PEC, parietal epithelial cell; CD-ICB, collecting duct type B intercalated cell; PST, proximal straight tubule cell; DCT, distal convoluted tubule cell; PODO, podocyte cell; ENDO, endothelium cell; PCT, proximal convoluted tubule cell; FIB, fibroblast cell; MES, mesangial cell; LEUK, leukocyte cell; KEGG, the Kyoto Encyclopaedia of Genes and Genomes; GO, gene ontology; UMAP: uniform manifold approximation and projection.

Functional enrichment of DEGs in scRNA-seq data

The results of the GO analysis indicated that the majority of up-regulated DEGs were enriched in the plasma membrane, playing a role in cell adhesion and extracellular matrix organization (Figure 1D and Table 1). In contrast, the down-regulated DEGs were mostly associated with urogenital system and kidney development in the apical regions of cells, as determined by GO enrichment analysis (Figure 1E and Table 1). Additionally, KEGG enrichment analysis revealed that the up-regulated genes were primarily enriched in pathways related to cell adhesion (Figure 1D and Table 1), while the down-regulated genes were mainly associated with the AGE-RAGE pathway and mineral absorption (Figure 1E and Table 1).

Identification of five different patterns of RCD gene markers and the assessment of RCD enrichment

We compiled a list of gene biomarkers associated with five different patterns of regulated cell death by utilizing the GO, KEGG, and FerrDb v2 databases, as presented in Table 2. Next, we utilized these RCD-related gene signatures to conduct gene set variation analysis on RNA-seq datasets, including single-cell and bulk data. In the bulk dataset, higher GSVA scores for necroptosis and apoptosis were observed in DN samples, while the score for ferroptosis was relatively low (Figure 2A). Interestingly, the scRNA-seq dataset also demonstrated higher necroptosis scores in DN samples. Moreover, we identified a higher pyroptosis score in DN samples based on the scRNA-seq data (Figure 2B).

Figure 2. The identification of necroptosis-related biomarkers. (A) GSVA score of five different types of programmed cell death (PCD) in the bulk RNA-seq data (GSE96804). (B) GSVA score of five different types of programmed cell death based on the scRNA-seq data (GSE131882). (C) Heatmap showing the correlation between PCD traits and gene modules. (D) Venn diagram demonstrating the necroptosis-related biomarkers by intersecting necroptotic gene modules arising from WGCNA analysis, DEGs in the bulk RNA-seq data and DEGs in the scRNA-seq data. Scatter plots showing Pearson’s correlation analysis between necroptotic GSVA scores and EGF (E), PAG1 (F) and ZFP36 (G) respectively. (**P<0.01; ***P<0.001). Abbreviations: DN: diabetic nephropathy; GSVA: gene set variation analysis; WGCNA: weighted gene correlation network analysis; DEGs: differentially expression genes.

Identification of necroptosis-related gene modules

We employed the WGCNA algorithm on the bulk RNA-seq dataset to investigate the gene modules associated with necroptosis. As a result, we identified three significant gene modules, indicated as the yellow, blue, and cyan modules in Figure 2C. By combining these modules, we identified 5090 genes that are associated with necroptosis (Supplementary Table 3).

Identification of necroptosis-related biomarkers

We aimed to identify necroptosis-related gene markers in our study. Firstly, we generated Venn diagrams for the scRNA-seq DEGs, the DEGs obtained from the bulk data (Supplementary Table 4), and the necroptosis-related gene modules (Figure 2D). By overlapping these gene lists, we identified three necroptosis-related gene markers: EGF, PAG1, and ZFP36. Next, we assessed the correlation between the expression levels of these biomarkers and necroptosis scores using Pearson’s correlation test. We found that EGF was negatively correlated with the necroptosis score (Figure 2E), whereas PAG1 was positively associated with the necroptosis score (Figure 2F). However, no statistically significant association was found between the expression levels of ZFP36 and the necroptosis score (Figure 2G).

Biomarker expression profiles

In the scRNA-seq dataset, EGF was mainly expressed at higher levels in distal convoluted tubule (DCT) cells of healthy subjects (Figure 3A, 3B). PAG1 was predominantly expressed in parietal epithelial cells (PEC) of patients with DN (Figure 3C, 3D). The expression levels of ZFP36 were significantly higher in cells of the proximal straight tubule (PST), fibroblasts (FIB), loop of Henle (LOH), and leukocytes (LEUK) of healthy subjects (Figure 3E, 3F).

Figure 3. Expression profiles of necroptosis-related biomarkers in the scRNA-seq data. A dimension reduction plot (A) and boxplot (B) showing the profiles of EGF expression in different cell types. The expression levels of PAG1 in scRNA-seq data as shown as a dimension reduction plot (C) and boxplot (D). ZFP36 expression profile in different types of cells in the GSE131882 as shown in a dimension reduction plot (E) and boxplot (F). (Wilcoxon’s test; *p<0.05; **p<0.01; ***p<0.001). Abbreviations: LOH, loop of Henle cell; CD-ICA, collecting duct type A intercalated cell; CT, connecting tubule cell; PEC, parietal epithelial cell; CD-ICB, collecting duct type B intercalated cell; PST, proximal straight tubule cell; DCT, distal convoluted tubule cell; PODO, podocyte cell; ENDO, endothelium cell; PCT, proximal convoluted tubule cell; FIB, fibroblast cell; MES, mesangial cell; LEUK, leukocyte cell; DN: diabetic nephropathy; UMAP: uniform manifold approximation and projection.

Moreover, based on the analysis of the bulk RNA-seq datasets, EGF and ZFP36 were expressed at significantly higher levels in the control group (p-value<0.001), whereas PAG1 was expressed at higher levels in subjects with DN (p-value<0.001) (Figure 4A). The ROC analysis results revealed that EGF, PAG1, and ZFP36 had high accuracy in discriminating between normal patients and those with DN, with AUC values of 0.877, 0.920, and 0.974, respectively (Figure 4B).

Figure 4. Expression profile of necroptosis-related gene markers in the bulk RNA-seq data (GSE96804) and validation in the test cohort (GSE142025). (A) Violin plot showing the expression levels of genes related to necroptosis in the bulk RNA-seq data. (B) Receiver operating characteristic curves of EGF, PAG1 and ZFP36 in the bulk RNA-seq data. (C) The expression levels of necroptosis-related biomarkers as validated by the external independent dataset and shown as a violin plot. (D) The performance of EGF, PAG1 and ZFP36 in the external independent dataset as determined by receiver operating characteristic curves. (**p<0.01; ***p<0.001, compared DN and control group).

Validation in the external independent cohort

In the external test cohort (GSE142025), we observed that EGF and ZFP36 were significantly overexpressed in normal subjects (p-value<0.001), while PAG1 was expressed at high levels in DN patients (p-value<0.01) with expression levels similar to those in the training cohort (Figure 4C). ROC analysis demonstrated that EGF and ZFP36 had strong discriminatory abilities for distinguishing between normal controls and DN patients (AUC = 0.865 for EGF; AUC = 0.977 for ZFP36), while PAG1 had a moderate discriminatory ability (AUC = 0.787) (Figure 4D).

Clinical features associations with the necroptosis-related gene markers

In this study, we investigated the potential of EGF, ZFP36, and PAG1 as biomarkers for chronic kidney disease (CKD). We found that EGF and ZFP36 were positively correlated with the estimated rate of glomerular filtration (P=2.3e-06 and R=0.83 for EGF; P=0.0082 and R=0.55 for ZFP36) (Figure 5A, 5G), indicating their potential to estimate the extent and progression of functional renal loss. However, PAG1 expression levels had no significant relationship with the glomerular filtration rate (P=1.8e-06 and R=-0.57 for PAG1) (Figure 5D).

Figure 5. Pearson correlation analysis of necroptosis-related biomarkers and clinical features. The scatter plots revealed a positive correlation between the expression level of the EGF gene (A) in DN patients and the glomerular filtration rate (GFR), while exhibiting a negative correlation with proteinuria (B) and serum creatinine (C). Additionally, the expression level of PAG1 was negatively correlated with GFR (D) and proteinuria (E), but positively correlated with serum creatinine (F). ZFP36, on the other hand, exhibited a positive correlation with GFR (G) but a negative correlation with proteinuria (H) and serum creatinine (I) levels.

Proteinuria, a known indicator of CKD, was negatively correlated with all three biomarkers (R=-0.56 for EGF; R=-0.8 for PAG1; R=-0.96 for ZFP36, all p-value <0.05) (Figure 5B, 5E, 5H). Furthermore, we found that the expression levels of EGF and ZFP36 were negatively correlated with serum creatinine levels in DN patients (P=0.0051 and R=-0.65 for EGF; P=0.0013 and R=-0.66 for ZFP36) (Figure 5C, 5I), whereas the levels of blood creatinine showed an opposite trend with regards to the expression of PAG1 in patients with DN (P=5.6e-06 and R=0.54 for PAG1) (Figure 5F).

Activation of necroptosis in HK-2 cells with high-glucose treatment

The necroptotic activation of HK-2 cells under high-glucose conditions was evaluated by western blotting and RT-qPCR analysis of RIP1, RIP3, and MLKL expression at the protein and mRNA levels. Our findings indicated that high glucose exposure led to a significant upregulation of RIP1, RIP3, and MLKL at both the protein (Figure 6A, 6B) and mRNA (Figure 6C) levels.

Figure 6. Activation of necroptotic and necroptosis-related biomarkers in HK-2 cells stimulated by high glucose. (A) Western blot assays were used to detect the expression levels of necroptotic proteins in HK-2 cells under different treatments. (B) Bar charts show the grey values obtained from the western blots shown in (A). (C) Reverse transcription quantitative PCR analysis was used to quantify the expression levels of necroptotic gene markers in HK-2 cells treated with NG, NG+M and HG, respectively. The protein and mRNA expression levels of EGF, PAG1 and ZFP36 were detected by western blot analysis (D, E) and reverse transcription quantitative PCR (F) in HK-cells when treated with NG, NG+M and HG. (#P < 0.05, ##P < 0.01, compared with NG group). Abbreviations: NG: normal glucose; NG+M: NG + Hyperosmotic medium; HG: high glucose; DN: diabetic nephropathy.

Expression of necroptosis-related biomarkers in HG-stimulated HK-2 cells

To evaluate the expression of EGF, ZFP36, and PAG1, we employed western blotting and RT-qPCR methods (Figure 6D–6F). We found that the expression levels of EGF were significantly lower in the HG group at both the mRNA and protein level when compared to the normal glucose treatment. Conversely, the expression levels of PAG1 were significantly higher in the HG group than in the NG group. However, there was no statistically significant difference in the expression levels of ZFP26 between the NG and HG groups.

Data collection and processing

To identify relevant datasets for our study, we used the keywords “Homo sapiens,” “diabetic kidney disease,” and “diabetic nephropathy” to search the Gene Expression Omnibus (GEO) database. We then carefully selected only those datasets that exclusively included patients with nephropathy caused by type 2 diabetes. Ultimately, we identified three datasets that met our criteria. The first was a single-cell RNA sequencing dataset (GSE131882), which included kidney tissue samples from three patients with DN and three healthy controls [39]. The second was a bulk RNA-seq dataset (GSE96804) that contained kidney tissue samples obtained from 40 patients with DN and 21 healthy controls [40]. Finally, we used the GSE142025 dataset as an external validation cohort, which consisted of 28 samples from patients with DN and 9 samples from healthy controls [41]. In Supplementary Table 5, we provide detailed clinical profiles of all the patients included in the three datasets.

In this study, we employed the Seurat package to analyze the single-cell sequencing data, as per the established guidelines [42]. We applied the following quality control standards to ensure the reliability of our analysis: (1) exclusion of cells with a gene count of less than 200 or more than 7500, (2) removal of cells with more than 25% of mitochondrial genes, and (3) elimination of double cells. The integrated samples were subjected to the Harmony algorithm and principal component analysis for linear downscaling and cell clustering. To visualize the data, we applied the non-linear downscaling method, “uniform manifold approximation and projection” (UMAP) algorithm. The function “FindAllMarkers” was utilized to identify specific genes in each cluster. Further, cell-specific markers mentioned in previous literature were used to annotate the clusters.

To analyze the bulk-RNA seq data, we first aligned the reads to the hg38 reference genome, and then transformed the raw read counts data to transcripts per million (TPM) for further analysis. Differential gene expression analysis was performed using the “limma” package [43]. We considered genes with a |log2-fold change (FC)| ≥ 1 and Benjamini and Hochberg adjusted p-value < 0.05 as differentially expressed genes (DEGs). These DEGs were subsequently used for downstream analyses.

DEGs in scRNA-Seq and functional enrichment

In accordance with the cluster annotation results, we used the ‘FindMarkers’ function to investigate the DEGs between DN and control individuals in various clusters. Genes with a p-value < 0.05 were identified as DEGs. Afterward, we employed the online enrichment analysis tool DAVID (http://david.ncifcrf.gov) and Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) databases to identify the biological functions and pathways of the DEGs [44]. We regarded the enrichment terms with a Benjamini and Hochberg adjusted p-value of less than 0.05 as statistically significant.

Identification of regulated cell death biomarkers and gene set variation analysis (GSVA)

To investigate the mechanisms underlying regulated cell death (RCD) in diabetic nephropathy (DN) samples, we searched the KEGG, GO and FerrDb v2 databases to identify RCD markers of five types, including autophagy, necroptosis, pyroptosis, apoptosis and ferroptosis [45]. To assess the enrichment of these RCD gene sets in the RNA sequencing datasets, we used gene set variation analysis (GSVA) to generate GSVA scores for each gene set, which evaluate the change in activity of the associated gene set [46]. We then created RCD gene sets from the collected RCD-related key genes and used them to assess the enrichment of various types of RCD. After calculating the GSVA scores of different RCD patterns, we selected RCD patterns with statistical differences (p-value<0.05) in RNA-seq data (both single-cell and bulk) for further investigation.

Identification of key RCD-related gene modules

We employed the weighted gene co-expression network analysis (WGCNA) algorithm, a data reduction and unsupervised classification method, to build a gene co-expression network and confirm co-expressed genes and modules associated with RCD traits [47]. Gene modules were identified using hierarchical clustering trees, and topological overlap matrix-based hierarchical clustering was used to detect them. Pearson correlation coefficients were calculated to determine the correlation of each module with RCD features, and modules significantly associated with RCD traits (p-value<0.05) were identified. The genes within these modules were then exported for further analysis.

Identification of key RCD-related biomarkers

In order to identify RCD-related biomarkers in DN samples, we first intersected the DEGs identified in both the single-cell and bulk RNA-seq data with the gene modules acquired by WGCNA. Subsequently, we investigated the expression profiles of these biomarkers in different types of renal cells using the single-cell sequencing dataset. Additionally, we also investigated the expression profiles of these candidate gene markers in bulk RNA-seq data. To assess the diagnostic ability of these gene markers, we employed the “ROCR” package to generate receiver operating characteristic (ROC) curves.

Validation in the external test cohort

To further evaluate the robustness of the biomarkers obtained from the training cohort, we assessed their expression levels in different groups and generated receiver operating characteristic (ROC) curves based on an independent external dataset (GSE142025).

Association between gene markers and clinical features

The Nephroseq database (http://www.nephroseq.org) is a valuable resource for investigating the association between key genes and important clinical renal parameters such as serum levels of creatine (SCr), glomerular filtration rate (GFR), and proteinuria (PRO) in DN patients. Therefore, we utilized the Nephroseq database to assess the correlation between key genes and these clinical parameters. Pearson correlation analysis was performed to investigate the relationship between key genes and GFR, SCr, and PRO in DN patients.

Cell culture and treatment

HK-2 cells were obtained from the National Collection of Authenticated Cell Cultures and cultured in DMEM (HyClone; Cytiva, USA) supplemented with 10% FBS (HyClone; Cytiva, USA) at 37° C and 5% CO2. When cells reached 70%-80% confluency, the medium was replaced with serum-free culture medium and incubated for 24h. The cells were then divided into three groups: high glucose (HG, 30 mmol/L glucose), normal glucose (NG, 5.5 mmol/L glucose), and normal glucose with hyperosmotic medium (NG+M, 5.5 mmol/L glucose + 24.5 mmol/L mannitol). These groups were cultured for 48h, and the proteins were extracted for subsequent experiments. The cell assays were performed in triplicate.

Western blot analysis

Protein was extracted from HK-2 cells using RIPA lysate (cat. no. P0013 Beyotime, China), and the protein concentration was quantified using the BCA kit (cat. no. P0010; Beyotime, China). A total of 50μg of protein was mixed with loading buffer in a 1:5 ratio, denatured at 100° C for 10 min, and then separated by SDS-PAGE and transferred onto a PVDF membrane (Bio-Rad Laboratories, Inc., USA) at 90V for 2h. The membranes were blocked using 5% skimmed milk powder for 2h at 37° C and then incubated overnight at 4° C with various primary antibodies, including anti-MLKL (cat. no. ab184718, Abcam, USA), anti-RIP3 (cat. no. ab209384, Abcam, USA), anti-RIP1 (cat. no. ab202985, Abcam, USA), anti-EGF (cat. no. ab218831, Abcam, USA), anti-ZFP36 (cat. no. orb39206, Biorbyt, UK), anti-PAG1 (cat. no. orb67086, Biorbyt, UK), and anti-beta actin (cat. no. ab115777, Abcam, USA), all diluted by 1:1000. The membranes were washed three times with TBST for 10 min each and then incubated with HRP-labelled secondary antibodies (cat.no. ab97051, Abcam, USA) diluted by 1:1000 for 1h. Finally, the ECL method was used for visualization, and the optical density of the bands was quantified using ImageJ (version 1.8). The experiments were performed in triplicate.

Reverse transcription quantitative PCR (RT-qPCR) analysis

Lastly, we evaluated the mRNA expression levels of RIP1, RIP3, MLKL, EGF and ZFP36 in HK-2 cells using specific primers (see Supplementary Table 6). Total RNA was extracted from HK-2 cells using TRIzol reagent (Invitrogen, USA), and cDNA was synthesized using PrimeScript RT Master Mix (Takara, Japan). Real-time PCR was performed using SYBR Premix Ex Taq II (Takara, Japan) under the following conditions: (1) initial denaturation at 95° C for 30 s; (2) annealing at 95° C for 55 s and 60° C for 40 s, for 40 cycles; and (3) melting at 95° C for 15 s, 60° C for 60 s, and 95° C for 15 s. The endogenous control was β-actin, and data were analyzed using the 2-∆∆Cq method.

Statistical analysis

In this study, all statistical analyses were performed using R software (version 4.2). For microarray data, we applied a screening threshold of |log2-fold change (FC)| ≥ 1 and a Benjamini and Hochberg adjusted p-value < 0.05. For single cell data, we considered a p-value < 0.05 and an absolute value of log2FC greater than 0.25 as screening conditions. The GSVA scores were compared between DN and control groups for each RCD gene set using appropriate statistical tests like Wilcoxon rank sum test. Only those RCD types showing significantly different GSVA scores between groups (p<0.05) were considered enriched and taken for further analysis. In WGCNA, only gene modules showing significant correlation with regulated cell death (p<0.05) were considered associated with the phenotype. To compare the expression of hub biomarkers between diabetic nephropathy and control groups, we used the Wilcoxon test.

Data availability statement

All public datasets used in this study can be downloaded in GEO by searching the accession ID provided in this article. The data that support this finding are available within the article or its supplementary materials.