Identification and validation of the biomarkers related to ferroptosis in calcium oxalate nephrolithiasis

Data analysis

The GSE73680 dataset was employed in this study, consisting of 6 normal renal papillary tissues from individuals without a history of kidney stones (referred to as the “Normal” group) and 29 RPs from patients with calcium stones (referred to as the “Plaque” group). Then, 484 ferroptosis-related genes were achieved from the FerrDb database (http://www.zhounan.org/ferrdb/current/) (Supplementary Table 1). A Venn plot was generated to display the intersection between these ferroptosis-related genes that were differentially expressed in the GSE73680 dataset. Among these, we identified 61 genes as differentially expressed ferroptosis-related genes (DEFERGs) based on the criteria of |Log2 FoldChange|≥1 and P < 0.05 (Supplementary Table 2).

Weighted gene co-expression network analysis

We first computed pairwise correlations between all genes in the dataset to establish the basis for constructing the co-expression network. Then, using the formula amn = |cmn|β, where cmn represents Pearson’s correlation between gene m and gene n, we constructed a weighted adjacency matrix. The parameter β determined the soft threshold power value. Next, we calculated different β values to assess the scale-free topology fit index and chose the one that yielded the best scale-free network. Following, the topological overlap measure (TOM) matrix quantifies the network connectivity of a gene by summing its adjacency values with all other genes in the network. Subsequently, genes with similar expression profiles were grouped into gene modules with average linkage hierarchical clustering. Afterwards, the gene significance (GS), and module membership (MM) were calculated for each gene with clinical traits. The GS represents the correlation between gene expression and the clinical trait, while the MM quantifies the degree of connectivity of a gene within its module. The nephrolithiasis-related modules were determined by assessing the coefficient of determination (R²) and statistical significance (P < 0.05). Finally, we identified 22 key DEFERGs within the nephrolithiasis-related modules.

Machine learning

To screen the hub DEFERGs, a machine learning algorithm was employed. Specifically, the least absolute shrinkage and selection operator (LASSO) regression approach was utilized with “glmnet” R package. Through the LASSO analysis, two DEFERGs were identified as the hub DEFERGs based on their importance and contribution to the model. To assess the discriminatory capacity of the two hub DEFERGs in distinguishing between Plaque group and Normal group, a receiver operating characteristic (ROC) curve was built to discriminate between these groups.

Consensus clustering and functional enrichment analysis

RPs samples were grouped under an unsupervised hierarchical clustering analysis with the “ConsensusClusterPlus” R package. To ascertain the optimal number of clusters, consensus matrix plots, cumulative distribution function (CDF) plots and the relative change in area under the CDF curve were employed. Additionally, principal component analysis was used to evaluate the disparities between clusters. The “limma” R package was utilized to identify differentially expressed genes (DEGs) between clusters, with a threshold of |log2(fold change)| > 3 and False discovery rate (FDR) < 0.01. Gene Ontology (GO) enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis were conducted with an adjusted P < 0.05. To further explore significant functional differences between clusters, we employed gene set enrichment analysis (GSEA) and Gene Set Variation Analysis (GSVA) with “clusterProfiler” R package by FDR < 0.25 and P < 0.05, allowing for the evaluation of enrichment of predefined gene sets or pathways within the clusters.

CaOx stone rat model and treatment

Male Sprague-Dawley rats aged 5-6 weeks were divided equally with n = 6 into three groups: the control group, the CaOx group, and the ferrostatin-1 (Fer-1) group. The control group had access to water ad libitum, while the CaOx group received 1% ethylene glycol (EG) in their drinking water. The Fer-1 group was treated with 1% EG in their drinking water and received an additional gavage administration of 0.25 mg/kg/d of Fer-1. Fer-1 was initially dissolved in dimethyl sulfoxide (DMSO) and then further diluted with 0.9% NaCl solution. The resulting solution had a final concentration of 0.25 mg/ml for Fer-1 and 0.1% for DMSO, which was administered via gavage to the Fer-1 group. Meanwhile, the control group and the CaOx group were gavaged with equal volumes of 0.9% NaCl solution containing 0.1% DMSO. After four weeks, blood plasma and renal tissues were collected for further analysis.

Histologic analysis

Various histologic staining techniques were employed to visualize specific components and proteins in the renal tissues. Hematoxylin and eosin (HE) staining was performed to observe the general structures of the renal tissues. Von Kossa (VK) staining was used to detect calcium deposits, while diaminobenzidine (DAB)-enhanced Prussian blue (PB) staining was utilized to identify iron deposits. Trivalent iron in tissues can form a blue precipitate with potassium ferricyanide, indicating a high content of iron elements in the tissue. Subsequent addition of DAB triggers an oxidation-reduction reaction between DAB and the precipitated trivalent iron, forming a brown compound, which improves the specificity of detection. Masson staining was performed to visualize collagen fibers and assess fibrosis in the renal tissues. Immunohistochemistry (IHC) and immunofluorescence (IF) were applied to detect the expression levels and localization of specific proteins in the renal tissues. Detailed information about the antibodies can be found in Supplementary Table 3. A blind histological evaluation of all sections stained was undertaken by two independent investigators to prevent bias. Either the person conducting the histologic analysis or the individual providing the samples is unaware of certain information.

Cell culture, cell counting kit-8 (CCK-8) and lactate dehydrogenase (LDH) assay

The HK-2 cell line, sourced from the Chinese Academy of Sciences in Shanghai, China, was subjected to various treatments. All the initial HK-2 cells were seeded equally in each case after counting with trypan blue. Since the effect of 1 mM oxalate on HK-2 cells for 24h was the most significant, we chose 1 mM concentration to carry out the follow-up experiments [18, 19], while a subset of HK-2 cells, designated as the Fer-1 group, received a 2-hour pretreatment with Fer-1 (2 μM). After accurately calculating the number of cells with trypan blue, 3 × 10³ HK-2 cells from different groups were seeded in a 96-well plate. Subsequently, 10 μl of CCK-8 solution was added at the designated time and read at 450 nm absorbance. Cellular LDH activities were evaluated with the LDH Assay Kit (Cat# C0016, Beyotime, Shanghai, China) as per the manufacturer’s guidelines.

Measurement of LPO and ROS levels

Malondialdehyde (MDA), glutathione (GSH), and superoxide dismutase (SOD) activities were quantified by MDA kits (Cat# A003-1), GSH kits (Cat# A005), and SOD kits (Cat# A001-1) (JianCheng, Nanjing, China), respectively. The levels of LPO of HK-2 cells in different groups were evaluated by Lipid Peroxidation Probe BDP 581/591 C11 (Cat# L267, Dojindo Laboratories, Japan) according to the manufacturer’s instructions. The assessment of reactive oxygen species (ROS) production was conducted using DCFH-DA and DHE staining, in accordance with the guidelines provided by the manufacturer, and was subsequently measured through microfluorimetry detection.

Quantitative real-time PCR (qRT-PCR) and molecular docking

The total RNA from the renal tissues was extracted using the TRnaZol RNA Kit (Cat# M5102, NCM, Suzhou, China). qRT-PCR assays were performed using the SYBR Green PCR system (Cat# R311, Vazyme, Nanjing, China) independently replicated three times. The sequences of primers were also included in Supplementary Table 3. The structure of LAMP2 (PDB ID: 2MOF) and MDM4 (PDB ID: 3FDO) was obtained from the Protein Data Bank (PDB, https://www.rcsb.org/), while the structure of Fer-1 (PubChem CID: 4068248) was achieved from PubChem database (https://pubchem.ncbi.nlm.nih.gov/). Molecular docking was carried out using AutoDock Vina (version 1.1.2) for all calculations in this project as described in our previous study [18].

Statistical analysis

The statistical analysis was conducted using R software (version 4.1.2) and GraphPad Prism (version 8.0). The Wilcox test, Student’s t-test, or ANOVA test were employed to compare two or more independent test series, respectively.

Data availability statement

The original contributions presented in the study are included in the article/Supplementary Materials. Further inquiries can be directed to the corresponding author.

Identification of DEFERGs

We first used the Venn plot to screen 409 ferroptosis-related genes that have intersections with the GSE73680 dataset (Figure 1A). With the threshold of |Log2 FoldChange|≥1 and adjusted P < 0.05, there are 61 ferroptosis-related genes identified as DEFERGs, and 33 of them were significantly upregulated in Plaque group, while 28 of them were significantly downregulated (Figure 1B). The heatmap showed the expression of these 61 DEFERGs in GSE73680 dataset (Figure 1C). An exploration into the protein-level connections between DEFERGs was conducted via the PPI network (Figure 1D). The word-cloud map showed the scale of the adjacent nodes of DEFERGs, reflecting the importance of PTEN, IFNG, SQSTM1, CAV1, MAPK1, PIK3CA, ADIPOQ, and MDM4 (Figure 1E). The adjacent nodes of DEFERGs were listed in Supplementary Table 4.

Figure 1. Identification of DEFERGs. (A) The Venn plot shows the number of ferroptosis-related genes that have intersections with the GSE73680 dataset (n=409). (B) The volcano plot of all DEFERGs between Normal and Plaque group (|Log₂ FoldChange|≥1, adjusted P < 0.05). (C) The heatmap of DEFERGs between Normal and CaOx group. (D) PPI networks of the 61 DEFERGs. (E) The word-cloud map shows the scale of the adjacent nodes of DEFERGs, reflecting the importance of genes.

Screening of the DEFERGs correlated with nephrolithiasis via WGCNA

We employed WGCNA to identify modules associated with nephrolithiasis. We set the soft threshold power value β = 8 to construct a scale-free network (Figure 2A). Subsequently, a hierarchical clustering tree was generated and modules were detected by dynamic tree cutting. Four color modules were identified (Figure 2B), among which the turquoise module demonstrated a significantly positive correlation (r = 0.31, P = 0.03) with nephrolithiasis (Figure 2C). Gene expression matrix of the turquoise module via WGCNA analysis was listed in Supplementary Table 5. Further, 22 DEFERGs in the turquoise module were identified with the Venn plot, which might be related to nephrolithiasis (Figure 2D).

Figure 2. Screening the hub DEFERGs via WGCNA and machine learning. (A) The power value β, set at 8, is indicated by the red line, ensuring a robust R² > 0.9. (B) The clustering dendrogram illustrates the amalgamation of gene co-expression modules, with each color signifying a distinct module. (C) The heatmap displays the correlation between modules and clinical traits, featuring correlation coefficients and corresponding P-values at the intersections of rows and columns. (D) The Venn plot illustrates the overlapping genes between DEFERGs and those found in the nephrolithiasis-related module (turquoise). (E, F) LASSO regression of the 22 intersected DEFERGs. (G) Coefficients between the two hub DEFERGs and the clinical portraits.

Screening of the hub DEFERGs via machine learning

We further employed LASSO regression analysis to identify potential hub genes associated with CaOx nephrolithiasis using the expression profiles of 22 DEFERGs (Figure 2E, 2F). LASSO regression analysis revealed two variables, LAMP2 and MDM4, as diagnostic markers for CaOx nephrolithiasis (Figure 2G). Subsequently, we observed a significant decrease in the expressions of LAMP2 and MDM4 in Plaque group (Figure 3A). To assess the diagnostic performance of LAMP2 and MDM4, we constructed ROC curves and obtained area under the curves (AUC) values of 0.793 and 0.805, respectively (Figure 3B, 3C), indicating that LAMP2 and MDM4 gene signatures have potential diagnostic values.

Figure 3. Identification of two clusters based on hub DEFERGs. (A) The bean plot shows the differential expression of hub DEFERGs between the normal tissues and RPs. (B, C) ROC curve of hub DEFERGs in RPs’ diagnosis. (D, E) The consensus cumulative distribution function (CDF) plot (D) and the relative change in area under the CDF curve (E). (F) The consensus matrix plot is displayed for k = 2 clusters. (G) The PCA plot visually represents the distribution of these two clusters. (H) Boxplots display the expression patterns of hub DEFERGs in cluster 1 and cluster 2. (I) GSVA analysis of biological pathways is conducted between the two distinct clusters, with red indicating activated pathways and blue indicating inhibited pathways.

Identification of two clusters based on the hub DEFERGs

The consensus clustering analysis of the expression profiles of LAMP2 and MDM4 revealed two distinct subtypes of CaOx nephrolithiasis. The stability of the clustering results was confirmed by the CDF curve (Figure 3D) and the CDF Delta area curve (Figure 3E). The consensus matrix plot showed that Plaque group could be divided into two clusters, namely cluster 1 (n = 10) and cluster 2 (n = 19) (Figure 3F). The detailed clusters based on the hub DEFERGs were shown in Supplementary Table 6. The PCA plot further demonstrated a clear distinction between the two clusters (Figure 3G). Additionally, the boxplot indicated that cluster 1 had higher expression of LAMP2 and MDM4 than cluster 2 (Figure 3H). As shown in Figure 3I, GSVA results revealed that cluster 1 was enriched in lipid metabolism, amino acid metabolism, bile acid metabolism, kidney metabolism, signal transduction, cell junctions, and immune response, which exert a major influence on CaOx nephrolithiasis.

Functional distinctions between two clusters

A total of 4277 differentially expressed genes (DEGs) were identified between clusters (Supplementary Table 7). Compared to cluster 1, 4067 DEGs were upregulated and 210 DEGs were downregulated in cluster 2 (Figure 4A). G protein-coupled receptor signaling pathway, transmembrane signaling receptor activity, and humoral immune response were enriched in GO terms (Figure 4B). Lysosome, Focal adhesion, and Chemical carcinogenesis-reactive oxygen species were enriched in KEGG pathways (Figure 4C). As shown in Figure 4D, the GO analysis results, positive regulation of ion transport, regulation of cytosolic calcium ion concentration, cellular calcium ion homeostasis, calcium ion homeostasis, and cellular divalent inorganic cation homeostasis were significantly enriched in biological processes (BP); keratin filament, intermediate filament, intermediate filament cytoskeleton, ion channel complex, and mitochondrial respiratory chain complex I were significant involved in cellular components (CC); ion channel activity, cation channel activity, transmitter-gated ion channel activity, ligand-gated ion channel activity, and metal ion transmembrane transporter activity were significant included in molecular functions (MF). KEGG analysis demonstrated that the DEGs were enriched in cAMP signaling pathway, and Maturity onset diabetes of the young. These results suggest that the DEGs mainly participate in ion transport such as iron and calcium, which is closely related to ferroptosis and CaOx nephrolithiasis.

Figure 4. Differentially expressed genes (DEGs) screening and pathway enrichment analysis between the two clusters. (A) The volcano plot of all DEGs between two clusters. (B, C) GSEA analysis reveals GO terms (B) and KEGG terms (C) enriched by DEGs between the two clusters. (D) Enriched GO terms and KEGG pathways of DEGs. (E) Expression of the markers of ferroptosis in the two clusters. (F) Correlation of LAMP2 and MDM4 with the markers of ferroptosis.

Differences of ferroptosis markers between clusters

We further investigated the relationship between ferroptosis markers in two clusters, including SLC7A11, CHAC1, SLC40A1, TF, TFRC, FTH1, GPX4, HSPB1, and NFE2L2 (Figure 4E). The expression of ferroptosis suppressors (SLC7A11, GPX4, HSPB1, and NFE2L2) was higher in cluster 1 compared to cluster 2, while the expression of ferroptosis driver (CHAC1) was lower in cluster 1. Consistent with these, the heatmap revealed the positive correlation of LAMP2 and MDM4 with the ferroptosis suppressors, as well as the negative correlation of LAMP2 and MDM4 with the ferroptosis driver (Figure 4F). These results indicate that LAMP2 and MDM4 may inhibit ferroptosis in CaOx nephrolithiasis.

Validation of the hub DEFERGs in CaOx nephrolithiasis rat models

To investigate the expression of LAMP2 and MDM4 in CaOx nephrolithiasis, the EG-induced CaOx nephrolithiasis rat models were constructed. The renal function results of rats showed that serum levels of blood urea nitrogen (BUN), serum creatinine (Scr), uric acid (UA), and Fe³⁺ in the CaOx model group were significantly elevated (Table 1). HE, VK, and Masson staining demonstrated the presence of tubular dilation and injury, CaOx crystals, and fibrosis in the renal tubules of the CaOx rats, which was reversed by the administration of Fer-1 (Figure 5A). PB staining showed significant iron accumulation in the CaOx group, and IHC results displayed a decrease in GPX4 expression in the CaOx group (Figure 5A). As depicted in Figure 5B, the Fer-1 group and the control group exhibited notably elevated levels of SOD and GSH in contrast to the CaOx group, whereas the MDA and Fe²⁺ levels displayed an inverse trend.

Figure 5. Fer-1 reduces CaOx crystal deposition and kidney injury in CaOx nephrolithiasis rat models. (A) Pathological sections including HE, VK, and Masson staining show the degree of kidney injury, crystal deposition, and kidney fibrosis; PB staining shows the Fe³⁺, while IHC staining shows GPX4 protein expression in different groups (magnification×20; scale bar, 200 μm). (B) SOD, MDA, GSH, and Fe²⁺ levels are detected in different groups (mean ± SD, n = 6, *P < 0.05, **P < 0.01).

TUNEL assay revealed that the apoptotic nuclear changes were markedly greater in the CaOx group than in the control and the Fer-1 group (Figure 6A). DHE fluorescence staining displayed the highest concentrations of intracellular oxidant species in the CaOx group (Figure 6A). Double-labelling IF results demonstrated that the protein expression of LAMP2 and MDM4 was attenuated in the CaOx group, which was in agreement with the bioinformatic analysis (Figure 6B). These findings demonstrate that the CaOx nephrolithiasis rat models were established successfully and Fer-1 treatment significantly reduced kidney injury, iron accumulation, and LPO in CaOx nephrolithiasis, suggesting the involvement of ferroptosis, as well as LAMP2 and MDM4 in CaOx nephrolithiasis.

Figure 6. Fer-1 reduces CaOx-induced LPO and apoptosis, as well as LAMP2 and MDM4 expression in CaOx nephrolithiasis rat models. (A) TUNEL and DHE staining techniques are employed to evaluate apoptotic damage and the production of reactive oxygen species (ROS) in renal tissues (magnification×40; scale bar, 400 μm). (B) Representative immunofluorescence images show LAMP2 (green) and MDM4 (red) detection in renal tissues (magnification×40; scale bar, 400 μm).

Validation of the hub DEFERGs in oxalate-induced HK-2 cell models

The cytoprotective effect of Fer-1 against cell injury induced by oxalate was evaluated using CCK-8 and LDH assays (Figure 7A). Furthermore, qRT-PCR (Figure 7B) and western blot (Figure 7C) analysis demonstrated that the mRNA and protein expressions of LAMP2 and MDM4 were decreased in the oxalate-treated group compared to the control group but were reversed after treatment with Fer-1. The binding energy results of molecular docking of Fer-1 with LAMP2 and MDM4 were -82.44 kcal/mol (LAMP2, Figure 7D), and -5.526 kcal/mol (MDM4, Figure 7E). Next, the DCFH-DA and DHE fluorescence intensity was observed to increase in the oxalate-treated group, but decreased following Fer-1 treatment (Figure 7F).

Figure 7. Fer-1 suppresses oxalate-induced cytotoxicity and LPO, as well as LAMP2 and MDM4 expression in oxalate-induced HK-2 cells. (A) Cell viability and LDH release after different treatment (mean ± SD, n = 3, *P < 0.05, ***P < 0.001). (B) qRT-PCR assays show GPX4, LAMP2, and MDM4 expression in HK-2 cells with different treatment (mean ± SD, n = 3, *P < 0.05, **P < 0.01, ***P < 0.001, compared with Oxalate group). (C) The LAMP2 and MDM4 protein expression of HK-2 cells in different groups. (D, E) Schematic diagram of Fer-1 docking with LAMP2 (D) and MDM4 (E). (F) Reactive oxygen species (ROS) generation in HK-2 cells was assessed using DHE and ROS staining (magnification×40; scale bar, 100 μm).

Consistently, the results of LPO probe revealed that HK-2 cells treated with oxalate exhibited significantly elevated LPO level, but was also reduced after Fer-1 treatment (Figure 8A). Similarly, the MDA level was elevated in the oxalate-treated group, but the SOD and GSH levels were reduced. These findings indicate the potential of Fer-1 to protect human renal epithelial cells against oxalate-induced ferroptosis injury (Figure 8B). These data were consistent with bioinformatics screening, strongly suggesting that LAMP2 and MDM4 were potential biomarkers in CaOx nephrolithiasis.

Figure 8. Fer-1 suppresses oxalate-induced cytotoxicity and LPO in HK-2 cells. (A) Measuring cellular LPO with the C11 BODIPY 581/591 fluorescent probe in HK-2 cells with different treatment (magnification×40; scale bar, 100 μm). (B) SOD, MDA, and GSH activity are detected in HK-2 cells with different treatment (mean ± SD, n = 3, *P < 0.05, **P < 0.01, ***P < 0.001, compared with Oxalate group).