This article is mentioned in:

Introduction

Diabetic nephropathy (DN) represents 30–50% of chronic kidney disease cases, affecting 285 million individuals globally (1). DN is one of the most severe microvascular complications of diabetes (2) and a leading cause of end-stage renal disease in both developing and developed countries, creating social and economic burdens (3,4). The key pathological characteristics of DN include tubulointerstitial fibrosis, the thickening of the glomerular basement membrane, accumulation of the mesangial matrix and nodular glomerulosclerosis (5), accompanied by proteinuria, hypertension and a progressive decline in renal function. Microalbuminuria is a key indicator for the evaluation of the progression of DN in clinical practice (6); however, not all patients with DN exhibit notable proteinuria; thus, it may be unreliable to assess prognosis or disease severity based solely on proteinuria levels (7). Treatment strategies aim to attenuate renal damage progression by controlling blood glucose levels or intraglomerular pressure (8,9). Hyperglycemia is a key hallmark of DN progression, and intensive glycemic control can delay its progression (10). Glucagon-like peptide-1 receptor agonists, sodium-glucose co-transporter-2 inhibitors, dipeptidyl peptidase-4 inhibitors and thiazolidinediones decrease hemoglobin A1c levels and the risk of developing microvascular complications (11,12). Hypertension is a key risk factor for the onset and development of DN. Angiotensin-converting enzyme inhibitors and angiotensin II receptor blockers are the preferred first-line anti-hypertensive treatments for patients with both diabetes and hypertension (13,14). Due to the heterogeneity among patients with DN, not all individuals benefit from these treatment approaches (15). Consequently, it is important to identify novel diagnostic and therapeutic targets to control the onset and progression of DN.

Epithelial-mesenchymal transition (EMT) is the process by which epithelial cells transform into mesenchymal cells, and is key in development and wound healing, while pathologically contributing to fibrosis and cancer progression (16). EMT is classified into three types: Type I is associated with embryonic and organ development, type II is associated with inflammation and fibrosis and type III is involved in cancer progression (17). The assessment of EMT includes EMT-associated transcription factors, such as Snail, Twist and zinc-finger E-box binding homeobox 1, as well as biomarkers, including α-smooth muscle actin (α-SMA), vimentin, epithelial cadherin (E-cadherin) and tight junction protein-1. Previous studies have demonstrated that in a high-glucose environment, human renal proximal tubular epithelial cells (HK-2) and renal tissue of patients with DN exhibit decreased E-cadherin mRNA and protein expression, along with increased expression of Snail, vimentin and α-SMA at both the mRNA and protein level, confirming the occurrence of the EMT process in the kidneys of patients with DN (18,19). In addition, studies have revealed fibrotic changes in the tubulointerstitium and glomeruli of patients with DN (20,21). Fibrosis is characterized by the accumulation of myofibroblasts (22), where tubular epithelial cells and podocytes undergo EMT to become myofibroblasts, leading to glomerulosclerosis and interstitial fibrosis (23). The inhibition of EMT in epithelial cells prevents chronic kidney injury and fibrosis (20), therefore, understanding the role of EMT-associated genes in the progression of DN to renal fibrosis is key.

The aim of the present study were to identify key EMT-related genes involved in the pathogenesis of DN to provide new insights for the diagnosis and treatment of DN. In the present study, based on the GSE96804 training set, key EMT-related genes were identified by determining differentially expressed genes (DEGs) and performing weighted gene co-expression network analysis (WGCNA). A protein-protein interaction (PPI) network of key genes was constructed and hub genes were identified using the Degree, Maximum Neighborhood Component (MNC) and Maximal Clique Centrality (MCC) algorithms in Cytoscape, then identifying optimal diagnostic genes through machine learning [least absolute shrinkage and selection operator (LASSO) regression analysis and support vector machine recursive feature elimination (SVM-RFE)]. The expression levels of diagnostic genes were validated in the validation sets GSE47183 and GSE104948, and receiver operating characteristic (ROC) curves were then generated to examine the diagnostic efficacy. To investigate the distribution of diagnostic genes in cell subpopulations, the Kidney Integrative Transcriptomics (K.I.T.) database was used for the single-cell analysis. Using the CIBERSORT algorithm, immune cell infiltrates in DN samples were calculated and the association between diagnosis genes and immune cells was assessed. Finally, reverse transcription-quantitative PCR (RT-qPCR) was used to confirm the expression of diagnosis genes.

Materials and methods

Human renal biopsy samples

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.

Animal experiments

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 CO2 vol/min was used for euthanasia. Mice were monitored for signs of respiration and pupillary reflex for 5 min following CO2 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).

Cell culture

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% CO2. For FN1 gene silencing, cells were seeded at 0.9×105/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.

Data sources

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).

Identification of DEGs

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.

WGCNA and selection of module genes

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 R2 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.

Connectivity map (CMap) 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.

PPI network construction

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.

Clustering analysis to identify the potential subtypes of DN

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).

Machine learning

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.

Hub gene validation and ROC curve analysis of GEO datasets

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).

Single-gene Gene Set Enrichment Analysis (GSEA)

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.

Single nucleus RNA sequencing

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.

Immune infiltration analysis

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.

RT-qPCR

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′.

Statistical analysis

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).

Results

Identification of DEGs

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).

Figure 1. Identification of DEGs in dataset no. GSE96804 from the Gene Expression Omnibus database. (A) Volcano plot of DEGs in dataset GSE96804. (B) Heatmap of DEGs. (C) UMAP illustrating the different groups. DEG, differentially expressed gene; DN, diabetic nephropathy; UMAP, uniform manifold approximation and projection.

WGCNA

WGCNA was performed to identify modules most strongly associated with DN. β=11 (scale-free R2=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).

Figure 2. Construction of the co-expression network. (A) β=11 was selected as the soft threshold with the combined analysis of scale independence and mean connectivity based on top 25% of genes with the highest variance. (B) Gene hierarchy tree-clustering diagram. Colors represent different module genes. (C) Heatmap of the correlation between module genes and phenotypes. (D) Venn diagram of the module genes from WGCNA, differentially expressed in GSE96804 and EMT-related genes. EMT, epithelial-mesenchymal transition; WGCNA, weighted gene co-expression network analysis.

Identification of candidate small molecular compounds for the treatment of DN

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).

Figure 3. Screening of potential small molecule compounds for the treatment of diabetic nephropathy via CMap database. (A) Top 10 compounds with the most significantly negative connectivity scores across 9 different cell lines. (B) Chemical structures of compounds. CMap, connectivity map; pc, principal component.

PPI network construction and important module analysis

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).

Figure 4. Screening hub genes using the STRING database by intersecting the top genes obtained from degree, MCC and MNC algorithms. (A) Protein-protein interaction network. Top 10 genes based on (B) degree, (C) MCC and (D) MNC. (E) Venn diagram of the genes screened from degree, MCC and MNC. MNC, maximum neighborhood component; MCC, maximal clique centrality.

Identification of potential subtypes of DN

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).

Figure 5. Potential molecular subtypes of diabetic nephropathy based on hub genes from dataset no. GSE96804 of the Gene Expression Omnibus. (A) Consensus matrix heatmap when k=2. (B) Representative CDF curves. (C) Score of consensus clustering. (D) Principal component analysis plot of the distribution of two clusters. CDF, cumulative distribution function; Dim, dimension

Identification of core genes through machine learning and validation

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).

Figure 6. Identification of diagnostic markers using machine learning to conduct validation. (A) Least absolute shrinkage and selection operator regression algorithm. (B) Support vector machine recursive feature elimination algorithm. (C) Expression of diagnostic markers in the GSE47183 and GSE104948 datasets from Gene Expression Omnibus. *P<0.05, **P<0.01. (D) ROC curve analysis of diagnostic markers in the GSE47183 and GSE104948 datasets. DN, diabetic neuropathy; AUC, area under the curve; CI, confidence interval; CV, cross validation; TPR, true positive rate; FPR, false positive rate.

GSEA

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).

Figure 7. Gene set enrichment analysis of FN1. Co-enriched (A) Gene Ontology terms and (B) Kyoto Encyclopedia of Genes and Genomes pathways.

Single nucleus RNA sequencing

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).

Figure 8. Single nucleus RNA sequencing from the Kidney Integrative Transcriptomics database. (A) Annotation of cluster subtypes. (B) Expression distribution of FN1. (C) Percent expression of FN1 in cell clusters. PCT, proximal convoluted tubule; CD, collecting duct; ICA, Type A intercalated cell; PEC, parietal epithelial cell; PC, principal cell; DCT, distal convoluted tubule; CT, connecting tubule; LOH, loop of Henle; PODO, podocyte; ENDO, endothelium; MES, mesangial cell; LEUK, leukocyte; tSNE, t-distributed stochastic neighbor embedding; FN1, fibronectin 1; CTRL, control.

Immune infiltration analysis

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).

Figure 9. Immune cell infiltration analysis based on dataset GSE96804 from the Gene Expression Omnibus database. (A) Infiltrating immune cells in DN and normal samples. (B) Spearman's correlation analysis of FN1 and immune cells. *P<0.05, **P<0.01 and ***P<0.001. DN, diabetic neuropathy; NK, natural killer; FN1, fibronectin 1; Cor, correlation.

Table I. Correlation analysis of fibronectin 1 and immune cells.

Table I.

Correlation analysis of fibronectin 1 and immune cells.

Immune cells	Correlation	P-value
B cells naive	−0.152730779	2.41×10−01
B cells memory	−0.599558608	3.29×10−07
Plasma cells	0.643217266	2.27×10−08
T cells CD8	0.065211452	6.18×10−01
T cells CD4 naive	0.553938893	3.64×10−06
T cells CD4 memory resting	0.298162117	1.96×10−02
T cells CD4 memory activated	−0.15596555	2.30×10−01
T cells follicular helper	−0.005988965	9.64×10−01
T cells regulatory	−0.048345671	7.11×10−01
T cells γ Δ	−0.077522704	5.53×10−01
NK cells resting	−0.271832858	3.41×10−02
NK cells activated	−0.377754714	2.69×10−03
Monocytes	−0.200384935	1.22×10−01
Macrophages M0	−0.12041248	3.54×10−01
Macrophages M1	0.608780767	1.94×10−07
Macrophages M2	0.415557753	8.67×10−04
Dendritic cells resting	0.766225281	6.21×10−13
Mast cells resting	−0.566155473	3.06×10−06
Mast cells activated	0.153979471	2.36×10−01
Eosinophils	−0.236738132	6.62×10−02
Neutrophils	−0.205340588	1.12×10−01

[i] NK, natural killer.

Validation of the expression of FN1 by RT-qPCR

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.

Figure 10. Validation of the expression of FN1 in the kidney of mice with DN and in patients with DN. Relative expression of FN1 mRNA in (A) kidney of mice (n=6) and (B) human kidney renal biopsy samples. Relative expression levels of (C) FN1, (D) E-cadherin and (E) vimentin in HK-2 cells (n=3). Data are expressed as the mean ± SD. *P<0.05. FN1, fibronectin 1; siRNA, small interfering RNA; DN, diabetic neuropathy.

Discussion

DN is a severe chronic complication of diabetes and exhibits morphological changes, including mesangial cell hypertrophy and proliferation, renal epithelial cell EMT and podocyte apoptosis (23,38). EMT is a process by which epithelial cells acquire mesenchymal features. Under physiological conditions, epithelial cells are connected by tight and adherens junctions and desmosomes, thereby maintaining the stability of cell polarity (39). Under pathological conditions, the induction of EMT transcription factors disrupts intercellular junctions, leading to transformation of epithelial into motile mesenchymal cells, which contribute to fibrosis (40). During the progression of DN, tubular epithelial cells and podocytes undergo EMT transformation into myofibroblasts, expressing fibroblast markers. Myofibroblast proliferation leads to glomerulosclerosis and interstitial fibrosis, exacerbating DN (41,42). The common final pathway of chronic kidney diseases, including DN, is renal fibrosis, with EMT being an initiating factor (43).

Studies have demonstrated that the activation of EMT-associated signaling pathways leads to development of renal fibrosis (44,45), while the inhibition of EMT attenuates renal fibrosis (46,47). EMT carries out a key role in the onset and progression of DN. Identifying key EMT-related genes may enhance the current understanding of the pathogenesis of DN and provide novel insight for its treatment. The present study retrieved transcriptomic datasets from the GEO database and identified nine key EMT-related genes using WGCNA and PPI networks. Based on the expression of these nine hub genes, consensus clustering analysis was performed on DN samples to explore potential molecular subtypes and reveal the heterogeneity of DN. Machine learning identified FN1 as the optimal diagnostic gene for DN. FN1 was validated using additional datasets and ROC curves were plotted to assess the diagnostic value of FN1. FN1 exhibited elevated expression in DN and demonstrated proficient diagnostic performance, a finding corroborated by previous studies (48,49).

FN1 encodes fibronectin (48) which is present in the extracellular matrix and participates in normal cell adhesion, proliferation, tissue development and wound healing (50,51). The overexpression of FN1 activates the TGF-β/PI3K/Akt signaling pathway, promoting bone fracture healing (52). FN1 is associated with the progression of various types of cancer, including squamous cell (53) and nasopharyngeal carcinoma (54) and ovarian (55) and renal cancer (56). Elevated expression of FN1 in tumor cells is negatively associated with the prognosis of patients with colorectal cancer (57). In addition, studies have indicated the key role of FN1, leading to glomerulosclerosis and fibrosis in chronic kidney disease (58,59). Fibronectin is a marker of mesenchymal cells (60). In the present study, single-gene GSEA was performed for FN1. As demonstrated by GO, FN1 is implicated in EMT, positioning it as a biomarker for this process (61). Knockdown of FN1 in tumor cells leads to increased expression of E-cadherin, while the expression of N-cadherin and vimentin decreases, indicating that FN1 knockdown inhibits the EMT process (62). As shown by KEGG pathway analysis, FN1 was predominantly enriched in the TGF-β signaling pathway. TGF-β1 is situated within FN fibrils and triggers EMT, while promoting fibrosis in DN (63–67). In the present study, compared with untreated HK-2 cells, the expression of E-cadherin in TGF-β-treated HK-2 cells decreased, while the expression of vimentin increased. These results indicated that TGF-β can initiate EMT. Additionally, FN1 is involved in modulating TGF-β-induced EMT, serving a key role in progression of DN (48,68). This conclusion was also verified in the present study. Silencing of FN1 led to upregulation of E-cadherin and downregulation of vimentin in HK-2 cells treated with TGF-β, indicating that FN1 was involved in the regulation of TGF-β-induced EMT.

The progression of DN is influenced by immune infiltration (69,70). CIBERSORT, a deconvolution method, was previously used to analyze immune cell gene expression patterns and determine the cellular composition within the tissue (34). The proportions of 22 immune cell types in DN samples were previously estimated, revealing a significant increase in M1 macrophages and resting dendritic cells (71–73). Activated macrophages serve a key role in the development of DN, leading to renal tissue matrix proliferation, glomerulosclerosis and irreversible pathological changes (74,75). Macrophages are classified into M1 and M2 types (76), with M1 macrophages serving a key role in antigen presentation and immune-inflammatory actions, while M2 macrophages release anti-inflammatory cytokines (77). A previous study confirmed the upregulation of the M1 macrophage marker CD86 in DN samples (73). Dendritic cells are antigen-presenting cells that activate and regulate T cell immune responses, serving a key role in the inflammation process in DN (78). Correlation analysis revealed that FN1 expression was positively correlated with plasma and CD4 naive and memory resting T cells, M1 and M2 macrophages and resting dendritic cells, while it was negatively correlated with memory B, resting and activated NK and resting mast cells. Specific T cell subsets have been found to conduct a key role in the development of DN (79). CD4+ T cells have been implicated in fibrogenesis, promoting renal fibrosis through mechanisms that involve myofibroblast aggregation and matrix deposition (80). A previous study using a mouse model of renal fibrosis demonstrated the upregulation of major histocompatibility complex class II (MHC II) expression in kidney tissue (81). MHC II molecules are essential for immune regulation, serving a key role in antigen presentation and T cell development. They facilitate the presentation of antigens to CD4+ T cells, thereby initiating adaptive immune responses (80). Therefore, FN1 may contribute to the development of DN by upregulating multiple types of immune cells.

The present study had limitations. First, the present study primarily relied on bioinformatics analysis, and although FN1 was associated with immune infiltration, experimental validation of this association was not conducted. Future studies incorporating in vitro and in vivo experiments are warranted to confirm the functional role of FN1 in immune regulation and its potential mechanistic involvement in disease progression. Second, the effect of FN1 on EMT was investigated only at the cellular level. Future studies should focus on generating FN1 knockout mouse models to explore the impact of FN1 silencing on EMT markers at both the mRNA and protein levels. Third, immunohistochemical or immunofluorescence analyses were not performed to assess the localization and expression of FN1 or EMT markers in renal tissue from diabetic mice or human patients. Moreover, the expression of EMT markers was not validated at the protein level by western blotting. Fourth, although the present results suggested that FN1 may regulate EMT through the TGF-β signaling pathway, the precise molecular mechanisms remain unclear. Further studies are required to elucidate how FN1 regulates the TGF-β pathway and contributes to EMT in DN. Finally, due to the limited number of available clinical samples, a reliable correlation analysis between FN1 expression and estimated glomerular filtration rate could not be performed. This limitation should be addressed in future studies with larger sample sizes to assess the clinical relevance of FN1 in kidney function decline.

In conclusion, FN1 was identified as the optimal diagnostic marker for DN. FN1 may serve key roles in the initiation and progression of DN by participating in EMT and upregulating various types of immune cell.

Supplementary Material

Supporting Data

Supporting Data

Supporting Data

Acknowledgements

Not applicable.

Funding

The present study was supported by The National Nature Science Foundation of China (grant. no. 82170734) and Shandong Provincial Natural Science Foundation, China (grant. no. ZR202102240178).

Availability of data and materials

The data generated in the present study may be requested from the corresponding author.

Authors' contributions

XW designed and supervised the study. KL analyzed data. KL and HT wrote the manuscript and performed bioinformatics analysis. HT revised the manuscript. HT and XW confirm the authenticity of all the raw data. All authors have read and approved the final manuscript.

Ethics approval and consent to participate

All experimental protocols for animal studies were approved by the Institutional Animal Care and Use Committee of School of Basic Medical Sciences, Shandong University (Jinan, China), and conducted in accordance with the National Institutes of Health Guide for the Care and Use of Laboratory Animals (approval no. ECSBMSSDU2023-2-157). The investigations of clinical specimens were conducted in accordance with the principles of the Declaration of Helsinki and were 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.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

References