The present study was conducted in accordance with the principles of the Declaration of Helsinki and approved by the Research Ethics Committee of Shandong University (approval no. ECSBMSSDU2023-1-157). Written informed consent was obtained from all patients prior to tissue collection. Renal biopsy samples (four male and two female patients) were collected from May 2023 to January 2024 as part of routine clinical diagnostics at the Department of Pathology, Qilu Hospital, Shandong University, Jinan, China. Inclusion criteria were as follows: i) Patients who underwent renal biopsy as part of routine clinical diagnostics and were diagnosed as DN by renal biopsy analysis; (2) availability of sufficient renal tissue for histological and molecular analysis; and (3) informed consent for use of biopsy tissue in research. Exclusion criteria included: (1) patients with known systemic infections or malignancies at the time of biopsy; (2) inadequate biopsy material and iii) absence of consent for research use. Control samples were obtained from the healthy kidney poles of individuals without kidney disease (three each male and female) who underwent tumor nephrectomy. The age range was 45–68 years.

For the present study, db/db (n=6) and db/m (n=6) mice were obtained from Beijing Vital River Laboratory Animal Technology Co., Ltd. All mice (3–5 mice/cage; 6–8 weeks old; weight 20–35 grams) were housed under standard specific-pathogen-free laboratory conditions, including a controlled temperature of 22±2°C, relative humidity of 50–60%, and a 12-h light/dark cycle. Mice had free access to water and a standard chow diet (Beijing Keao Xieli Feedstuff Co., Ltd.). The water and cages were autoclaved and cages with corncob bedding were changed three times/week. As female mice are resistant to obesity and type 2 diabetes (24), only male mice were selected. Mice were randomly assigned to groups, and the investigators were blinded to group allocation during surgery and outcome evaluations. A volume of 30% vol/min CO₂ vol/min was used for euthanasia. Mice were monitored for signs of respiration and pupillary reflex for 5 min following CO₂ administration to confirm death. Death was also confirmed by observing the absence of heartbeat. The protocol complies with AVMA guidelines and institutional ethical standards (25).

Human tubule epithelial cells (HK-2 cells) were obtained from American Type Culture Collection and cultured in a 50:50 mixture of DMEM/F12 supplemented with 10% fetal bovine serum (both Gibco, Thermo Fisher Scientific, Inc.) under 37°C, 5% CO₂. For FN1 gene silencing, cells were seeded at 0.9×10⁵/well in 6-well plates and cultured in the DMEM/F12 without antibiotics. Upon reaching 30% confluency, cells were transfected with short interfering RNA (siRNA; cat. no. AM16708; Thermo Fisher Scientific, Inc.; assay ID, 10826) for FN1 or scramble control siRNA (cat. no. AM4611; Thermo Fisher Scientific, Inc). A total of 10 nmol/l siRNA was used in this experiment. siRNA was delivered using Lipofectamine 3000^® reagent (Thermo Fisher Scientific, Inc.) following the manufacturer's protocol. Cells were transfected at 37°C for 6 h, after which the medium was replaced with fresh culture medium. Twenty-four hours after transfection, cells were treated with recombinant human TGF-β1 (10 ng/ml; Sigma-Aldrich; Merck KGaA) for an additional 24 h at 37°C.

Gene Expression Omnibus (GEO) database (ncbi.nlm.nih.gov/geo) esd searched using the formula ‘[‘diabetic nephropathies’(MeSH Terms) OR diabetic kidney disease(All Fields)] AND ‘Homo sapiens’(porgn) AND [‘Expression profiling by array’(Filter) AND (‘2013/01/01’(PDAT): ‘2023/01/01’(PDAT)]’. Datasets GSE96804 (26), GSE47183 (27) and GSE104948 (28) were included in the present study. GSE96804 was the training set with GSE47183 and GSE104948 serving as the validation sets. GSE96804 was downloaded from the GPL17586 platform [HTA-2_0 Affymetrix Human Transcriptome Array 2.0 (transcript (gene) version)], comprising 20 normal and 41 glomerular samples from patients with DN. GSE47183 was downloaded from the GPL14663 platform [Affymetrix GeneChip Human Genome HG-U133A Custom CDF (Affy_HGU133A_CDF_ENTREZG_10)], including 14 glomerular samples from patients undergoing nephrectomy for tumors and seven from patients with DN. GSE104948 was downloaded from the GPL22945 platform (HG-U133_Plus_2 Affymetrix Human Genome U133 Plus 2.0 Array [CDF: Brainarray HGU133Plus2_Hs_ENTREZG_v19)], consisting of three glomerular samples from patients undergoing tumor nephrectomy and seven glomerular samples from patients with DN. EMT-associated genes were obtained from the dbEMT 2.0 database (dbemt.bioinfo-minzhao.org/download.cgi).

Uniform manifold approximation and projection (UMAP) analysis of GSE96804 was conducted using umap (version 0.2.10.0) and ggplot2 packages (version 3.5.1) in R software (version 4.3.2) to visualize sample grouping. The limma package (version 3.58.1) was used to analyze GSE96804, identifying DEGs with P<0.05 and |log fold-change|≥1. Volcano plots and heatmap of DEGs were generated using the ggplot2 and pheatmap packages (version 1.0.12), respectively.

Using the WGCNA package (version 1.73) in R, a gene co-expression network was constructed to investigate the association between genes and phenotypes. The top 25% of genes with the highest variance were selected for analysis. The soft threshold was set to a power of 11, and a scale-free R² of 0.9 was applied to the weighted adjacency matrix, which was then transformed into a topological overlap matrix (TOM). ‘TOMType’ was set to ‘unsigned’. Modules were identified using the TOM-based dissimilarity metric (1-TOM) with a minimum module size of 30, based on hierarchical clustering. Each module was assigned a color. Among the 15 gene modules, the blue, pink and yellow modules, whose correlation with DN was ≥0.4, were selected for further analysis.

CMap (clue.io) is a gene expression profile database that captures gene expression changes induced by various perturbations, enabling the exploration of associations between diseases, genes and small molecule compounds. Upregulated EMT-related genes were analyzed using the CMap database to identify potential small molecule drugs for the treatment of DN. Small molecules with negative connectivity scores suggest potential to reverse EMT gene expression in DN, highlighting their therapeutic value. The top 10 compounds with the most significant negative connectivity scores were identified as potential candidates.

To explore protein interactions, pathways and co-expression, the STRING database (cn.string-db.org/) and Cytoscape software (version 3.10.0) (cytoscape.org/) were used to construct a PPI network of key EMT-related genes. The cytoHubba plugin was applied to identify significantly interacting genes. Hub genes were determined as the intersection of the top 10 genes identified by the Degree, MCC and MNC algorithms.

Based on the expression profiles of the nine hub genes most strongly associated with DN, the ConsensusClusterPlus package (version 1.66.0) was used for consensus clustering to unveil distinct subtypes of DN. The clusterAlg was ‘hc’ and distance was ‘pearson’. Scatterplot3d (version 0.3.44) and factoextra packages (version 1.0.7) were employed for principal component analysis (PCA).

To filter the optimal diagnostic genes for DN, LASSO and SVM-RFE were used. LASSO logistic regression analysis is a regression methodology proficient at selecting variables to enhance predictive accuracy (29). It serves as both a technique for regression and variable selection, employing regularization to improve predictive accuracy and enhance the comprehensibility of statistical models (30). SVM-RFE analysis is a supervised machine learning technique that identifies optimum diagnostic genes by progressively discarding feature vectors produced by SVM (31). The glmnet (version 4.1.8) and e1071 packages (version 1.7.16) were used to conduct LASSO and SVM-RFE analyses, respectively. The intersection genes of LASSO and SVM-RFE were considered optimal genes in the diagnosis of DN.

GSE47183 and GSE104948 were used as validation datasets. Box plots were constructed using the ggplot2 package to assess the expression of core genes. The area under the ROC curve was employed to evaluate the diagnostic efficacy of core genes for DN, with the ROC curve generated using the pROC package (version 1.18.5).

GO analysis enables a more comprehensive understanding of gene product functional annotation, whereas KEGG pathway analysis not only facilitates annotation of gene functions but also elucidates the roles of genes within diverse signaling pathways (32). To investigate the potential functions of core genes, GSEA-GO and GSEA-KEGG were performed using the org.Hs.eg.db (version 3.18.0) and clusterProfiler packages (version 4.10.1). Genes in the GSE96804 dataset were ranked based on Spearman's correlation coefficients with the core gene. P<0.05 was applied to identify enriched pathways.

The K.I.T. database (Ben Humphreys' Laboratory at the University of Washington; humphreyslab.com/SingleCell/) (33) was used to analyze single nucleus RNA data of core genes within DN samples and visualize the outcomes. The dataset included a total of 12 types of cell: Proximal convoluted tubule; collecting duct; type A intercalated cell; parietal epithelial cell; principal cell; distal convoluted tubule; connecting tubule; loop of Henle; podocyte; endothelium; mesangial cell and leukocyte.

CIBERSORT is a deconvolution algorithm designed to measure immune cell infiltration (ICI). Infiltration of 22 cell types was measured within the gene expression profile of GSE96804 (34). ICI outcomes in the samples were visualized utilizing the ggplot2 package. Spearman's correlation analysis was performed to ascertain the correlation between the core gene and the infiltrating immune cells, utilizing the ggcorrplot and ggplot2 packages.

Total RNA was isolated from kidney tissue of mice or renal biopsy samples using TRIpure (BioTeke Corporation) according to the manufacturer's instructions. Using a PrimeScript™ RT reagent kit (Takara Bio, Inc.), cDNA was synthesized with 1 µg RNA as a template. The following temperature protocol was used for RT: 85°C for 5 sec for RT; 37°C for 15 min to inactivate the reverse transcriptase. RT-qPCR was performed using a SYBR-Green qPCR kit (Takara Bio, Inc.) according to the manufacturer's instructions. The thermocycling conditions were as follows: Initial denaturation at 95°C for 10 min, followed by 40 cycles of 95°C for 30 sec, 58°C for 30 sec and 72°C for 20 sec. The housekeeping gene β-actin served as an internal control to normalize RNA quantity and quality across samples. Gene expression fold changes, normalized to β-actin, were calculated using the 2^−ΔΔcq method (35). The following primers were used: Mouse FN1 forward, 5′-CCCTATCTCTGATACCGTTGTCC-3′ and reverse, 5′-TGCCGCAACTACTGTGATTCGG-3′; β-actin forward, 5′-GGCTGTATTCCCCTCCATCG-3′ and reverse, 5′-CCAGTTGGTAACAATGCCATGT-3′ and human FN1 forward, 5′-ACAACACCGAGGTGACTGAGAC-3′ and reverse, 5′-GGACACAACGATGCTTCCTGAG-3′; E-cadherin forward, 5′-GCCTCCTGAAAAGAGAGTGGAAG-3′ and reverse, 5′-TGGCAGTGTCTCTCCAAATCCG-3′; vimentin forward, 5′-AGGCAAAGCAGGAGTCCACTGA-3′ and reverse, 5′-ATCTGGCGTTCCAGGGACTCAT-3′ and β-actin forward, 5′-CTCACCATGGATGATGATATCGC-3′ and reverse, AGGAATCCTTCTGACCCATGC-3′.

All experiments were independently repeated ≥3 times. Data are expressed as the mean ± SEM or SD. Statistical analyses were performed using R software (version 3.6.1) and GraphPad Prism software (version 8.0; Dotmatics). Normality assumption of the data distribution was assessed using Kolmogorov-Smirnov test. Comparisons between two groups were conducted using unpaired two-tailed Student's t-test for normally distributed data and the Mann-Whitney rank sum test for non-normally distributed data. Differences between >2 groups with a single variable were analyzed using one-way ANOVA, followed by Tukey's post hoc test. No power analysis was performed to determine the sample size. The sample size was based on our previous studies (36,37).

After analyzing the GSE96804 dataset, 624 DEGs were identified, comprising 284 up- and 340 downregulated genes (Table SI; Fig. 1A and B). The UMAP plot revealed a distinct separation between the DN and normal samples (Fig. 1C).

WGCNA was performed to identify modules most strongly associated with DN. β=11 (scale-free R²=0.9) was selected based on scale independence and mean connectivity (Fig. 2A). A total of 15 co-expression gene modules (CGMs) were generated (Fig. 2B). The correlation between DN and GCMs was calculated (Fig. 2C), with the blue, pink and yellow modules selected for subsequent analyses. The intersection (25 genes) of module genes, DEGs and EMT-related genes were obtained (Fig. 2D).

To explore potential small molecule drugs with therapeutic effects for patients with DN, the CMAP database was used to predict compounds capable of reversing EMT-related gene expression changes in DN. The top 10 compounds with the lowest connectivity scores were VU-0415374-1, marbofloxacin, propantheline, guanfacine, evoxine, PCO-400, β-CCP, SB-431542, CP-724714 and OM-137. These compounds were considered potential therapeutic agents for DN (Fig. 3A and B).

A PPI network (Fig. 4A) was constructed for core gene selection, employing MCC, MNC and Degree algorithms to identify the hub genes (Fig. 4B-D). In total, nine overlapping genes (FN1, MMP2, CCL2, TGFB2, MMP7, KDR, LTBP1, VCAN and TNC) were designated as hub genes (Fig. 4E).

To identify potential subtypes of DN, the present study performed consensus clustering on the 41 DN samples; k=2 provided the most stable clustering result. This was determined by evaluating the cumulative distribution function (CDF) curve, where k=2 showed the smallest increase in the area under the CDF curve, indicating higher stability (Fig. 5A and B). Furthermore, when k=2, the consistency score of each subtype is >0.9 (Fig. 5C). Therefore, the 41 DN samples were divided into Cluster 1 (n=21) and 2 (n=20; Table SII). The assignment of each sample to a cluster was determined by the consensus clustering algorithm, which calculates the frequency with which each sample pairs with others across multiple resampling iterations. Finally, PCA analysis revealed distinct clustering of these groups (Fig. 5D).

To determine the optimal diagnostic marker for DN, machine learning approaches were applied. Using LASSO, two genes were identified from the set of nine hub genes (Fig. 6A) and three genes were identified by employing SVM-RFE (Fig. 6B). The intersection gene (FN1) identified by both machine learning approaches was deemed the optimal diagnostic marker for DN. The expression levels of FN1 were validated, revealing elevated FN1 expression in DN samples from the GSE47183 and GSE104948 datasets (Fig. 6C). To evaluate the diagnostic efficacy of FN1, ROC curves were generated using the GSE47183 and GSE104948 datasets. The AUC value for the ROC curve in dataset GSE47183 was 0.898 and 1.000 in dataset GSE104948, indicating that FN1 may have significant diagnostic value (Fig. 6D).

As per Gene Ontology (GO) analysis, FN1 was primarily involved in ‘epithelial to mesenchymal transition’, ‘positive regulation of epithelial to mesenchymal transition’ and ‘regulation of epithelial to mesenchymal transition’ (Fig. 7A). As per Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis, FN1 was primarily involved in ‘TGF-beta signaling pathway’ (Fig. 7B).

By performing single nucleus RNA sequencing, the distribution of FN1 was determined in 12 cell groups (Fig. 8A); FN1 was primarily distributed in mesangial cells and the endothelium (Fig. 8B and C).

Using the CIBERSORT algorithm, the proportions of ICI in DN and normal tissue were calculated (Table SIII). There was a significant increase in the infiltration of M1 macrophages, M2 macrophages and resting dendritic cells and decrease in the infiltration of memory B, resting mast cells and eosinophils in the DN samples compared with normal tissue (Fig. 9A). To determine whether FN1 serves a key role in the pathogenesis of DN by regulating immune cells, the correlation between FN1 and immune cells was analyzed (Fig. 9B). FN1 was positively correlated with plasma cells (r=0.64), CD4 naive (r=0.55) and memory resting T cells (r=0.30), M1 (r=0.61) and M2 macrophages (r=0.42) and resting dendritic cells (r=0.77); however, it was inversely correlated with memory B (r=−0.60), resting (r=−0.27) and activated NK (r=−0.38) and resting mast cells (r=−0.57; Table I).

Relative mRNA expression of FN1 was elevated in mice with DN and in renal biopsy tissue samples from patients with DN (Fig. 10A and B). FN1 siRNA transfection was successful (Fig. 10C). Silencing of FN1 partially recovered E-cadherin expression in TGF-β treated HK-2 cells compared to cells treated with TGF-β alone (Fig. 10D). Notably, within the FN1 siRNA group, E-cadherin expression in TGF-β treated cells was significantly lower than in untreated controls, indicating a partial recovery rather than full normalization. FN1 silencing significantly ameliorated vimentin upregulation (Fig. 10E). Similarly, within the FN1 siRNA group, vimentin expression in TGF-β treated cells remained higher than in untreated controls, indicating a partial recovery rather than full normalization.