The cuproptosis-related gene ITGB6 and LTBP1 may be associated with diabetic kidney disease progression and immune cell infiltration

Data sources

Training set: The DKD dataset GPL17586 (Liu et al., 2018) was downloaded from the GEO database (https://www.ncbi.nlm.nih.gov/geo/), which utilizes chip sequencing (platform: GPL17586). The GPL17586 dataset was selected as the training set due to its comparatively larger sample size (40 DKD vs. 21 controls), well-annotated clinical phenotypes, and relevance to glomerular pathology in DKD. Kidney glomerular samples from 40 DKD patients were used as the DKD group, while glomerular samples from unaffected portions of the kidneys of 21 tumor nephrectomy patients were used as the control group.

Validation set: The validation set, GPL571 (Park et al., 2011), profiled on the GPL571 platform (Affymetrix Human Genome U133A 2.0 Array), was chosen as an independent cohort to evaluate the generalizability of our findings across different technical platforms. The GPL571 dataset was downloaded (platform: GPL571), including kidney glomerular samples from 9 DKD patients and 26 non-DKD patients as the DKD and control groups, respectively. Additionally, 36 cuproptosis-related genes were obtained from the literature (Chen et al., 2023b).

Data preprocessing subsection: To address potential batch effects between the GPL17586 and GPL571 datasets, which were generated on different microarray platforms (GPL17586 and GPL571, respectively), we applied the ComBat algorithm from the “sva” R package (version 3.46.0) to harmonize the expression data. The ComBat adjustment was performed using an empirical Bayes framework, with the dataset specified as the batch covariate and the disease group (DKD vs. control) preserved as a biological variable. This approach removed technical artifacts while retaining biologically relevant signals.

Differential expression analysis

The R programming language (version 4.3.2) and the Limma package (version 3.56.0) were used to identify differentially expressed genes (DEGs) between DKD and control groups in the GPL17586 dataset. Genes with an absolute log₂fold change (|log₂FC|) > 1 and an adjusted p-value (Benjamini-Hochberg method) < 0.05 were considered statistically significant. Visualization was performed using the R package ggplot2 (version 3.4.4) for volcano plots and ComplexHeatmap (version 2.16.0) for heatmaps.

ssGSEA analysis

Single-sample gene set enrichment analysis (ssGSEA), an extension of the gene set enrichment analysis (GSEA) method, was used to calculate enrichment scores for each sample based on the cuproptosis gene set. Wilcoxon rank-sum tests were applied to compare ssGSEA scores between the DKD and control groups, assessing the overall upregulation or downregulation of cuproptosis-related gene expression.

WGCNA analysis

Weighted gene co-expression network analysis (WGCNA) was performed using the WGCNA package (version 1.72-5) on gene expression data from 61 samples (DKD and control groups). A soft-thresholding power (β) of 5 was selected to ensure a scale-free topology (scale-free R² > 0.85). Co-expression modules were identified using a dynamic tree-cutting algorithm (minModuleSize = 30, deepSplit = 2, mergeCutHeight = 0.25). Modules with a correlation P-value < 0.05 and |correlation coefficient| > 0.3 with the cuproptosis trait were considered significantly associated. Module membership (MM) > 0.8 and gene significance (GS) > 0.5 were used as cutoffs for identifying hub genes within key modules.

Machine learning for key gene selection

To identify candidate key genes with diagnostic potential, we employed two machine learning algorithms: eXtreme Gradient Boosting (XGBoost) and Least Absolute Shrinkage and Selection Operator (LASSO) regression.

For the XGBoost model, we used the “xgboost” package (version 1.7.5) in R. The hyperparameters were tuned via 10-fold cross-validation repeated 3 times to optimize generalization performance. The grid search space included: max_depth (3, 6), eta (0.01, 0.1), gamma (0, 0.2), subsample (0.7, 1), colsample_bytree (0.7, 1), and min_child_weight (1, 5). The final model was trained using the area under the receiver operating characteristics (ROC) curve (AUC) as the evaluation metric. Feature importance was assessed by the “gain” metric, which reflects the relative contribution of each gene to model accuracy.

For the LASSO regression, we utilized the “glmnet” package (version 4.1.7) in R. The optimal regularization parameter λ was selected through 10-fold cross-validation that minimized the binomial deviance. The value of λ corresponding to the minimum cross-validation error (λ.min) was chosen to identify non-zero coefficients and retain robust features. Genes with non-zero coefficients under the optimal λ were considered candidate diagnostic biomarkers. Both models were trained and validated using the GPL17586 and GPL571 datasets. The resulting key genes were then intersected from each algorithm to obtain a set of high-confidence diagnostic markers that reached a consensus.

Expression validation

The Wilcoxon test was used to analyze the expression differences of candidate key genes between the DKD and control samples (P < 0.05) in the GPL17586 training set and the GPL571 validation set. Boxplots were generated with ggplot2 to visualize the results. Genes showing significant differences and consistent trends (P < 0.05) in both datasets were considered key genes.

ROC curve analysis

ROC curve analysis was conducted to evaluate the diagnostic performance of candidate genes in distinguishing between DKD and control samples. Using the R package “pROC” (version 1.18.0), ROC curves were generated and the area under the curve (AUC) was calculated for each gene in both the training set (GPL17586) and the independent validation set (GPL571). Genes with an AUC value greater than 0.7 were considered to exhibit good diagnostic ability. The ROC curves visualize the trade-off between sensitivity and specificity across different threshold values, providing a quantitative measure of each gene’s potential as a diagnostic biomarker for DKD.

Nomogram construction

A nomogram was constructed based on a multivariate regression model integrating key genes with diagnostic value for DKD. The total score was calculated based on the weighted scores of each key gene, and the DKD probability was predicted based on the total score.

Enrichment analysis

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were performed for DE-hub genes using the R package “clusterProfiler.” GO analysis included biological processes (BP), cellular components (CC), and molecular functions (MF). KEGG pathway analysis was used to identify major metabolic and signaling pathways involving the differentially expressed genes.

PPI network analysis

A protein-protein interaction (PPI) network was constructed for DE-hub genes using the STRING database and visualized with Cytoscape software. Topological properties of the network (e.g., node degree and centrality) were analyzed to identify key proteins.

Immune cell infiltration analysis

To investigate the relationship between key genes and the immune microenvironment, ssGSEA was performed using the R package “GSVA” on the GPL17586 dataset to calculate immune cell enrichment scores. Wilcoxon tests were used to compare immune cell infiltration differences between the DKD and control groups, and Spearman correlation analysis was conducted to examine relationships between immune cells and diagnostic genes.

In vitro mouse podocyte culture and differential gene expression analysis

(1) Cell culture: The mouse renal podocyte cell line MPC5 was obtained from Procell Life Science & Technology Co., Ltd. (catalog number: CL-0855; RRID: CVCL_0496). Cells were cultured in RPMI-1640 medium (Gibco, Waltham, MA, USA) supplemented with 10% fetal bovine serum (FBS; Procell, Bethel, CT, USA) and 1% penicillin/streptomycin (Beyotime, Shanghai, China) at 37 °C in a humidified 5% CO₂ atmosphere. High-glucose treatment was initiated when the cell density reached 80–90%. All cell experiments were conducted in accordance with relevant institutional guidelines and did not require additional ethical approval as the study exclusively used a commercially available cell line.

(2) High-glucose treatment: To study the effects of high glucose on LTBP1 and ITGB6 expression, MPC5 cells were divided into four groups: one control group (11.1 mM D-glucose concentration (Scrimieri et al., 2023)) and three high-glucose treatment groups (30 mM D-glucose concentration (Colombaioni et al., 2022)) for 24, 48, and 72 h, respectively.

(3) RNA extraction and qRT-PCR: After high-glucose treatment, total RNA was extracted from MPC5 cells using TRIzol reagent (Invitrogen, Waltham, MA, USA) following the manufacturer’s instructions. RNA concentration and purity were measured using a Nanodrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). cDNA was synthesized using a reverse transcription kit (Biosharp, Hefei, China). qRT-PCR was performed to detect the relative mRNA expression levels of LTBP1 and ITGB6 using specific primers and ACTIN as the internal control. Reactions were conducted on an ABI 7500 Real-Time PCR System using SYBR Green PCR Master Mix (Yesen, China). Each experimental condition included three independent biological replicates, with each biological replicate measured in technical triplicate. The relative expression levels were calculated using the 2^−ΔΔCt method, primer sequences were shown in Table 1.

(4) Statistical analysis: All experimental data are presented as mean ± standard deviation. A two-way analysis of variance (two-way ANOVA) was performed to assess the main effects of treatment (Control vs. HG) and time (24, 48, 72 h), as well as their interaction. When a statistically significant interaction was observed, post hoc tests with Bonferroni correction were further conducted to compare groups at each time point. All statistical analyses were carried out using GraphPad Prism version 9.0, with P < 0.05 considered statistically significant.

Differential gene expression and WGCNA analysis in the training set (GPL17586)

Differential gene expression analysis

To identify differentially expressed genes (DEGs) between DKD and control samples in the GPL17586 dataset, the “limma” package in R was used with thresholds of |log2FC| > 1 and adjusted P-value (P.adj) < 0.05. A total of 336 DEGs were identified, including 107 upregulated and 259 downregulated genes. A volcano plot was generated using the “ggpubr” package to visualize these results (Fig. 1A). The top 10 genes ranked by absolute log2FC values were selected for visualization in a heatmap created with the “ComplexHeatmap” package (Fig. 1B).

Candidate gene screening

To identify differentially expressed genes associated with cuproptosis, we intersected the 366 differentially expressed genes with the key genes from the two co-expression modules, defining these intersections as candidate genes. Venn diagrams were generated using the R package “VennDiagram” to visualize the results. The brown module contained 124 differentially expressed genes (Fig. 2A), while the blue module contained 28 differentially expressed genes (Fig. 2B).

Expression validation

To validate the expression of the key genes identified by machine learning, we compared the differential expression of candidate genes between DKD samples and control samples in both the training set (GPL17586) and the validation set (GPL571) using the Wilcoxon test (P < 0.05). The results showed significant differences in the expression of the candidate genes between DKD and control samples. The analysis results were visualized using box plots created with the R package “ggplot2.” Genes that showed significant differential expression (P < 0.05) and consistent trends in both the training and validation sets were defined as key genes.

As shown in Fig. 3A, the differential expression analysis of the key genes in the brown module from the training set (GPL17586) revealed significant differences for all genes. In contrast, the analysis in the validation set (GPL571) (Fig. 3B) showed no significant differences, and no expression was observed for KIAA1191 and MORN2. Figure 3C presents the differential expression analysis of the blue module in the training set (GPL17586), where all genes showed significant differences. Figure 3D presents the results from the validation set (GPL571), where significant expression differences were observed for ITGB6 and LTBP1, both of which demonstrated consistent upregulation in both the training and validation sets. Together, these results confirm ITGB6 and LTBP1 as key genes associated with DKD.

ROC analysis

To validate the diagnostic performance of the key genes for DKD, we conducted ROC curve analysis using the R package “pROC” in both the training set (GPL17586) and the validation set (GPL571), and calculated the area under the curve (AUC). As shown in Figs. 4A–4D, the AUC values for ITGB6 and LTBP1 in both the training and validation sets were greater than 0.7, indicating that both genes exhibit good diagnostic performance for DKD.

To further assess the impact of the key genes on the clinical probability of DKD, we constructed a nomogram based on all samples from the training set (GPL17586) using the R package “rms.” First, scores were assigned based on the diagnostic performance of each key gene for DKD, and the total score (Total Point) was calculated. The probability of DKD occurence was then inferred from the total score. A higher total score corresponds to a higher probability of disease occurrence. As shown in Fig. 4E, the model constructed using ITGB6 and LTBP1 had a total score of 113, resulting in a DKD prediction probability of 0.959, thereby demonstrating strong diagnostic potential.

The predictive ability of the nomogram model was further evaluated using calibration curves, decision curves (DCA), and ROC analysis. The calibration curve was assessed using the Hosmer-Lemeshow (HL) test. A P-value greater than 0.05 and a slope of the calibration curve close to 1 indicate good model fit. As shown in Fig. 4F, the P-value was 0.583, and the slope of the calibration curve was close to 1, suggesting that the model fits well.

Figure 4G presents the DCA decision curve, where the diagnostic model curve (blue) is above the gray baseline, and the net benefit (Net Benefit) is greater than 0, confirming that the model performs well in predicting DKD.

Finally, we assessed the predictive ability of the nomogram model in the training set using ROC analysis with the R package “pROC.” As shown in Fig. 4H, the AUC value was 0.89, demonstrating that the nomogram model possesses strong predictive capability.

Immune infiltration analysis

To explore the relationship between the key genes ITGB6 and LTBP1 and the immune microenvironment, we used the R package “GSVA” to calculate the enrichment scores of 28 immune cell types in DKD samples and control samples using the ssGSEA algorithm, based on the training set (GPL17586). As shown in Fig. 5A, a stacked bar plot was created to illustrate the proportion of immune cell infiltration between the two groups, with different colors representing different cell types, and the y-axis indicating the relative proportion of immune cells for each sample.

We applied the Wilcoxon test to compare immune cell infiltration between the two groups. As shown in Fig. 5B, 17 immune cell types exhibited significant differences (P < 0.05) between the groups, including activated B cells, CD56bright natural killer cells, central memory CD4 T cells, effector memory CD4 T cells, etc. To further investigate the relationship between ITGB6 and LTBP1 and these differential immune cells, Spearman correlation analysis was conducted using the R package “psych.” The results (Fig. 5C) indicated that ITGB6 and LTBP1 were significantly positively correlated with Activated B cells and Central memory CD4 T cells, and significantly negatively correlated with Neutrophils.

Key gene and related protein analysis

The PANTHER classification system is a comprehensive and annotated gene family database that categorizes proteins and their corresponding genes. To further investigate the classification of the key genes at the protein level, we utilized this database to analyze the key genes and obtain classification information for ITGB6 and LTBP1. The results are shown in Table 2.

Expression of ITGB6 and LTBP1 in renal podocytes under high glucose culture

To investigate the impact of a high-glucose environment on the expression of LTBP1 and ITGB6 in MPC5 cells, we exposed MPC5 cells to a high-glucose condition of 30 mM for 24, 48, and 72 h. Subsequently, we measured the relative mRNA expression levels of LTBP1 and ITGB6 using qRT-PCR.

As shown in Fig. 6, the expression of both LTBP1 and ITGB6 significantly increased with prolonged high-glucose treatment. For LTBP1, the relative expression levels at 48 h were significantly higher than the control group, with a P-value of 0.0012. Similarly, ITGB6 exhibited a significant time-dependent upregulation, with relative expression levels at 24, 48, and 72 h showing P-values of 0.035, 0.0014, and 0.0034, all significantly higher than the control group.

Limitations

This study indicates that cuproptosis-related pathways may contribute to DKD pathogenesis, with significant upregulation of ITGB6 and LTBP1 observed in DKD and their potential involvement in modulating immune microenvironment. Although our study has yielded significant findings, there are still several aspects that warrant further improvement. First, the absence of a mannitol control group limits our ability to exclude osmotic effects. Second, classical cuproptosis signaling pathways—including FDX1, LIAS, and DLAT protein lipoylation, as well as mitochondrial membrane potential alterations—require further investigation to better substantiate the link between cuproptosis and DKD. Finally, additional studies employing in vitro and in vivo DKD models, combined with Western blotting, immunostaining, and functional inhibition assays targeting ITGB6 and LTBP1, are necessary to elucidate their roles in cuproptosis and inflammatory signaling under hyperglycemic conditions.