Comprehensive Analysis of Fatty Acid Metabolism in Diabetic Nephropathy from the Perspective of Immune Landscapes, Diagnosis and Precise Therapy

Introduction

With the increasing prevalence of diabetes mellitus, diabetic nephropathy (DN) represents the most common complication of diabetes.^1,2 It is becoming the leading cause of end-stage renal disease worldwide, bringing a substantial economic burden to the public.^3–5 Excessive mesangial matrix production, renal hypertrophy, and fibrosis are the typical pathological changes in DN.^6–8 The pathogenesis of DN is complicated. Various risk factors, eg, deposition of advanced glycation end products (AGEs), boosting of the renin-angiotensin-aldosterone system (RAAS), and oxidative stress, are reported to exert critical effects on the initiation and evolution of DN.^9–12 Yet, the underlying mechanisms of DN remain largely unrevealed. Over the past few decades, only a few medications, including RAS inhibitors, Finerenone, sodium-glucose co-transporter-2 (SGLT-2) inhibitors, and glucagon-like peptide 1 (GLP-1) agonists, have been available for reducing the incidence of DN in populations with diabetes by controlling the glucose level or arterial pressure.^4,13,14 Nonetheless, the therapeutic drug strategy directly targeting DN is still limited. In addition, due to the heterogeneity and complicated mechanisms of diabetes-related injuries in the kidney, the existing diagnosis and therapy of DN remain controversial,^15,16 which are more complicated than other types of kidney diseases.¹⁷ Therefore, further exploring the pathogenesis is of great significance, and identifying novel biomarkers for early diagnosis, preferable classification, and precise treatment of DN is urgently needed.

It is well acknowledged that the kidney is an active metabolic organ and requires vast energy to maintain the homeostasis of the microenvironment.^18–20 A growing number of studies have considered fatty acid (FA) as the primary energy source to support kidney functional needs.¹⁸ Dysfunction of FA metabolism could be observed in the experimental DN mouse models, as indicated by excessive deposition of lipid and decreased expression of carnitine palmitoyl transferase 1A (CPT1A), which is a rate-limiting enzyme in the FA metabolism of the kidney.^21,22 Recent studies have highlighted that impaired FA metabolism may play an essential role in the pathogenic procedures of multiple kidney injuries.^18,23 Overexpression of renal CPT1A markedly ameliorates kidney injuries and fibrosis in experimental animal models by modulating the oxidative process of FA metabolism.²⁴ Consistently, inhibition of FA oxidation significantly restores the damaged barrier function of renal vessels in DN mice.²⁵ Moreover, fatty acid-binding protein 4 (FABP4), an important fatty acid transporter protein, is involved in apoptosis of mesangial cells in DN via regulating endoplasmic reticulum (ER) stress.²⁶ Furthermore, previous studies have also shown that abnormal metabolism of multiple nutrients, including impaired FA metabolism, can alter the immune microenvironment via different pathways during the pathological processes.²⁷ Both innate and adaptive immune responses have been recently found to participate in the pathogenesis of DN.²⁸ However, the regulatory role of FA metabolism and its association with immune responses in DN have not been fully elucidated, and it remains undetermined whether some factors associated with FA metabolism could become promising biomarkers for DN diagnosis or classification.

In the present study, we further clarified the role of FA metabolism-related genes in DN by integrating multiply bioinformatics analysis. Patients with DN were firstly classified into two groups via clustering analysis. We then investigated the differences in functional enrichment and immune landscape as well as clinical outcomes between these clusters and also discovered candidate drugs for precise medication. Importantly, novel biomarkers (PRKAR2B and ANXA1) were identified and a diagnostic model for predicting DN risk was further established based on the dysregulated FA metabolism-related genes by two machine learning algorithms. Our study will provide new insights into the role of FA metabolism in the pathogenesis and the classification of DN, and the early diagnostic biomarkers and precise therapy for DN.

Methods

Data Processing

Microarray datasets of glomerulus tissues from patients with DN and control groups, including GSE96804, GSE30528, and GSE104948, were taken from the Gene Expression Omnibus (GEO) repository. The integrated DN dataset, which contained 57 DN and 54 control samples, was generated after the batch effect correction of these three datasets using the combat function of the “SVA” package²⁹ in R software.

Access of Fatty Acid Metabolism-Related Genes

Fatty acid (FA) metabolism-related genes were downloaded from the fatty acid metabolism process (ontology gene sets) in the Molecular Signatures Database (MsigDB) (https://www.gsea-msigdb.org/gsea/msigdb).

Identification of differentially expressed fatty acid metabolism-related genes

The “Limma” package³⁰ in R software was applied to the identification of all differentially expressed genes (DEGs) (DN vs Control) with the cut-off values of adjusted p-value <= 0.05 and |log2 Fold Change (FC)| ≥ 0.585. Venn plot was then used to define differentially expressed FA metabolism-related genes by intersection of FA metabolism-related genes downloaded from the database and all DEGs in DN.

Molecular Typing of Fatty Acid Metabolism-Related Genes

The “ConsensusClusterPlus” package³¹ in R software was used for clustering analysis to identify potential subtypes of the DN on the basis of differentially expressed FA metabolism-related genes. The optimal k-value was selected using the consensus matrix (CM) and the maximum cumulative distribution function (CDF) index. The DEGs between the two molecular subclusters were screened out by “Limma” package with the threshold of adjusted p-value <= 0.05, |log2 FC| ≥ 0.585.

Weighted Gene Co-Expression Network Analysis (WGCNA) and critical module genes determination

The “WGCNA” package³² in R software was carried out to determine key modules associated with different molecular subtypes of the DN samples. After obtaining the optimal soft threshold power, constructing a scale-free co-expression gene network, and calculating the correlation coefficient between modules and subtypes, the module-trait relationships were then built. Modules with highly positive correlations with different molecular subtypes were selected for further analysis.

Functional Enrichment Analysis

The roles of genes in different molecular subtypes were revealed by annotation of Gene Ontology (GO) and analysis of Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and p-value < 0.05 was regarded as significant.

Immune Landscape Analysis

The “CIBERSORT” package³³ in R software was employed to evaluate the abundances of the infiltrated immune cells from the DN dataset. The Wilcoxon test was used to determine the differences in the number of 22 types of immune cells between controls and DN samples, together with different molecular subtypes, respectively. The correlation among 22 types of immune cells, and the relationship between differentially expressed FA metabolism-related genes, and the infiltrated immune cells, were calculated using Spearman’s rank correlation coefficient. Adjusted p-value < 0.05 was determined to be significant.

Potential Small-Molecular Drugs Prediction

Connectivity map (CMAP)³⁴ (https://clue.io) serves as a repository for gene expression data on the alteration of gene expression signatures following treatment of various small-molecular compounds, providing information for predicting potentially therapeutic drugs. Up-regulated key genes in Cluster B, based on the molecular subtyping of FA metabolism-related genes, were imported into the CMAP database to predict the potential small-molecular drugs for future therapeutic application. The top 10 drugs with the highest negative enrichment scores were chosen as potential pharmacological agents.

Screening of Diagnostic Biomarkers and Establishment of DN Risk Prediction Nomogram

To screen the candidate diagnostic biomarkers, the Least Absolute Shrinkage and Selection Operator (LASSO) together with Support Vector Machine-Recursive Feature Elimination (SVM-RFE) algorithms were performed by the “glmnet”³⁵ and “e1071”³⁶ packages in R software, respectively. The genes included in both algorithms were defined as potential diagnostic biomarkers for further construction of DN risk prediction model using the “rms” package³⁷ in R software.

Evaluation for Diagnostic Efficacy of the Prediction Model

The calibration curves, decision curve analysis (DCA), and the receiver operating characteristic (ROC) curves were applied to evaluate the performance of the nomogram in diagnosing patients with DN. Moreover, the ROC curve was also conducted to assess the efficacy of the nomogram to predict the risk of an early DN stage.

Correlation Analysis of Clinical Features

The “Nephroseq v5” tool (https://nephroseq.org/) was utilized to evaluate the correlation between the expression of the candidate diagnostic biomarkers and clinical characteristics in patients with DN.

Establishment of an Experimental Animal Model

Six-week-old of male C57BL/6J mice were obtained from GemPharmatech Co. Ltd, China. After adaptive feeding for a week, mice were randomly divided into two groups and given access to a normal chow diet (NCD, 10 kcal% fat) or high-fat diet (HFD, 60 kcal% fat), respectively. After four weeks of NCD or HFD, the HFD mice were induced by streptozotocin (STZ; Sigma-Aldrich, USA) via intraperitoneal (i.p.) injection in a dose of 40 mg/kg body weight for seven consecutive days, whereas NCD mice were received the same volume of vehicle buffer through i.p. injection. After two weeks of the injection, the blood samples were obtained from the tail vein of mice and the random blood glucose levels were detected by an ACCU-CHEK^® Performa glucometer (Roche, Manheim, Germany).³⁸ The blood glucose levels more than 16.7 mmol/L were regarded as the successful establishment of an experimental diabetic animal model and the corresponding mice were finally used in this study. Mice were sacrificed after 11 weeks of NCD or HFD, and the serum and urine samples as well as kidney tissues were harvested to perform further experimental analysis. Animal experimental protocols gained approval from the Institutional Research Ethics Committee at the Sun Yat-sen University (SYSU-IACUC-2022-001575). They were carried out following the relevant guidelines and the Guide for the Care and Use of Laboratory Animals (NIH publications Nos. 80–23, revised 1996).

Renal Function and Histology

Serum creatinine (SCr) and urinary creatinine (UCr) were measured by using the Creatinine Assay kit (Nanjing Jiancheng Bioengineering Institute, Nanjing, China). Blood urea nitrogen (BUN) was detected by the Urea Nitrogen Content Assay Kit (Solarbio, Beijing, China). Urine albumin protein levels were determined using the Albumin assay kit (Nanjing Jiancheng Bioengineering Institute). Paraffin-embedded kidney tissues were sectioned (4 μm) and subjected to Periodic Acid-Schiff stain (PAS) and Masson’s trichrome staining. The areas of glomerular and mesangial matrix were quantified using ImageJ software version 1.53e (National Institutes of Health, Bethesda, USA), respectively, and the mesangial matrix index was obtained as the ratio of the mesangial matrix area to the glomerular area. The areas of collagen deposition were calculated using ImageJ software based on 10 randomly chosen high-power fields per kidney section.

Protein Extraction and Western Blotting (WB)

The kidney tissues were lysed in RIPA buffer (Beyotime, Shanghai, China) containing protease and phosphatase inhibitors. The proteins in the lysates were resolved via sodium dodecyl sulfate- polyacrylamide gel electrophoresis (SDS-PAGE), and transferred onto the polyvinylidene difluoride (PVDF) membranes.³⁹ The membranes were then blocked by 5% skim milk at room temperature (RT) for an hour, followed by incubation of anti-PRKAR2B (28351-1-AP; Proteintech, Wuhan, China; diluted 1:500), anti-ANXA1 (21990-1-AP; Proteintech; diluted 1:5000) or anti-GAPDH (60004-1-Ig; Proteintech; diluted 1:20000) antibodies at 4°C for 16 hours. On the second day, the membranes were incubated with the horseradish peroxidase (HRP)-conjugated secondary antibody (SA00001-1, SA00001-2; Proteintech; diluted 1:10000) at RT for an hour. Protein levels were assessed by WB Chemiluminescence Detection. GAPDH was used as the endogenous reference and the grey values of protein bands were quantified by ImageJ software.

Real-Time Quantitative Polymerase Chain Reaction (RT-qPCR)

RNA extraction was carried out by the Trizol reagent (ThermoFisher Scientific, Darmstadt, Germany) according to previous study,⁴⁰ and reverse transcription was conducted with the QuantiTect Reverse Transcription Kit (Ruizhen Bio, Guangzhou, China). RT-qPCR was performed via the SYBR Green PCR kit (Ruizhen Bio). The 2^−ΔΔCt approach was adopted to calculate the relative mRNA expression of PRKAR2B and ANXA1. Primers for β-actin, PRKAR2B and ANXA1 were as shown as followings, mouse PRKAR2B-F: 5’-CCAGTAAGGGTGTCAACTTCG-3’, mouse PRKAR2B-R: 5’-GGACTCTGCATCGTCTTCCTC-3’; mouse ANXA1-F: 5’-AAGCAGGCCCGTTTTCTTGAA-3’; mouse ANXA1-R: 5’-GCAACATCCGAGGATACATTGA-3’; mouse β-actin-F: 5’-GGCTGTATTCCCCTCCATCG-3’, mouse β-actin-R: 5’-CCAGTTGGTAACAATGCCATGT-3’.

Statistical Analysis

Data processing and analysis were conducted in R software (version 4.2.1). GraphPad Prism version 9.0.2 was performed for statistical analysis in RT-qPCR analysis, and an unpaired t-test was performed to calculate the differences between the two groups. Results were shown as mean ± standard deviation (SD), and p-value < 0.05 was regarded as significance.

Results

Differential Expression Analysis of FA Metabolism-Related Genes in Patients with Diabetic Nephropathy

The flow diagram of this present research process is performed in Figure 1. In this study, three microarray datasets were obtained from the GEO repository and integrated after eliminating the batch effects, finally consisting of 57 DN and 54 control glomerular samples (Additional file: Figure S1A and B). With the threshold of adjusted p-value < 0.05 and |log2 FC| > = 0.585, 610 differentially expressed genes (DEGs), including 219 up-regulated and 391 down-regulated genes, were identified and presented in the volcano plot (Additional file: Figure S1C) and the expression of partial DEGs was shown in the heatmap (Additional file: Figure S1D). Among them, a total of 39 FA metabolism-related genes (8 up-regulated and 31 down-regulated genes) were found to be dysregulated in patients with DN by overlapping genes in FA metabolic process with DN-associated DEGs (Additional file: Figure S1E and F).

Figure 1 The research flow chart of this study.

Identification of Molecular Subtypes Based on FA Metabolism-Related Genes

Based on the 39 dysregulated FA metabolism-related genes, patients with DN were divided into two subtypes (Cluster A, n = 31 and Cluster B, n = 26) using the “Consensus Cluster Plus” package, as k = 2 was found to be comparatively stable in classification (Additional file: Figure S2A–C). WGCNA was employed to identify the key gene modules in Cluster A and Cluster B with the soft-threshold power of β = 7 (R² = 0.929, Additional file: Figure S2D). In all, ten gene modules were obtained, and the cluster dendrogram of the modules was then created (Additional file: Figure S2E and F). Furthermore, the correlation coefficient of molecular subtypes with gene modules was calculated, and yellow together with green modules had high positive association with Cluster A, whereas the turquoise together with magenta modules exhibited high positive relative to Cluster B (Additional file: Figure S2G). Thus, these four modules were chosen as the key modules for further analysis. By functional enrichment analysis of the biological process (BP), the genes included in Cluster A-associated modules were mainly correlated to “cell cycle”, “response to endoplasmic reticulum stress” and “integrin-mediated signaling pathway” (Additional file: Figure S3A). For Cluster B-associated modules, the included genes were mainly enriched in “inflammatory response”, “immune response” and “innate immune response” (Additional file: Figure S3B). KEGG analysis revealed that the Cluster A-associated modules were principally related to “protein processing in endoplasmic reticulum”, while in the Cluster B-associated modules “cytokine-cytokine receptor interaction” was the major enriched pathway (Additional file: Figure S3C and D).

Determination of the Differences in Functional Enrichment and Immunological Characteristics Between the Clusters Based on FA Metabolism-Related Genes

In order to further explore the differences in the biologically relevant function between the molecular clusters based on FA metabolism-related genes, 545 DEGs were screened out according to the criteria of adjusted p-value < 0.05 and |log2 FC| > = 0.585. Volcano plot and heatmap were presented to show the significant difference in gene expression patterns between the two subtypes (Figure 2A and B). Subsequently, functional enrichment analysis was employed to clearly understand the role of DEGs between these two FA metabolism-related groupings. Gene Ontology (GO) term analysis showed that the DEGs were functionally annotated. Of note, the DEGs were mainly enriched in “inflammatory response”, “immune response” and “extracellular matrix organization” upon biological process analysis (Figure 2C). Interestingly, the expression levels of pro-fibrotic genes (FN1, COL1A1, and ACTA2) were significantly higher in Cluster B than those in Cluster A. Meanwhile, compared to Cluster B, pro-inflammatory genes, such as CCL18, CCL2, and TNFRSF12A, showed markedly lower expression levels in Cluster A, suggesting the distinct status of the patients with DN (Figure 2D).

Figure 2 Functional enrichment analysis and immune cell infiltration analysis between these two subtypes. (A) The volcano plot of the DEGs between these two subtypes. (B) The heatmap of the DEGs between these two subtypes. (C) GO analysis of the DEGs upon the biological process. (D) The differential expression of FN1, COL1A1, ACTA2, CCL18, CCL2 and TNFRSF12A between these two subtypes. (E) The stacked histogram for the immune cell abundances between Cluster A and Cluster B of DN patients. (F) The violin diagram exhibiting the differences of immune cells between Cluster A and Cluster B. Red asterisks showing the increased number of immune cells infiltration in Cluster B, whereas blue asterisks indicating the reduction in Cluster B. *p-value < 0.05; **p-value < 0.01; ***p-value < 0.001; ****p-value < 0.0001. Abbreviation: ns, not significant.

Given that DEGs were closely associated with inflammatory and immune processes, the CIBERSORT algorithm was performed to clarify the characteristics of immune cell infiltration according to gene expression in immune cells between these two groupings. The proportions for 22 types of immune cells in each sample are presented in Figure 2E. Notably, five subpopulations of immune cells exhibited significant differences between cluster A and cluster B (Figure 2F). Samples in Cluster A had higher proportions of “Dendritic cells resting”, “Plasma cells”, and “T cells regulatory (Tregs)” than those in Cluster B. However, the lower abundances of “Macrophages M2” and “T cells gamma delta” were observed in Cluster A compared to Cluster B (Figure 2F).

Relationship Between the Differentially Expressed FA Metabolism-Related Genes and Immune Cell Infiltration

The immune landscape in kidney samples of patients with DN was dramatically changed (Figure 3A). The correlation analysis for 22 types of immune cells in DN samples showed that the proportion of “T cells CD8” was significantly negatively correlated to that of “T cells CD4 memory resting” (r = −0.73), and the abundance of “Mast cells activated” was negatively associated with that of “Mast cells resting” (r = −0.86) (Figure 3B). We further explored the relationship between the expression of differentially expressed FA metabolism-related genes and the abundance of immune cells (Figure 3C). Interestingly, the results indicated that “T cells gamma delta” was positively associated with PLA2G4A, SOX9, CD36, and C3, whereas “Macrophages M2” was positively correlated to APOC1 but negatively correlated to PTGDS (Figure 3D and E). In addition, infiltrated “Neutrophils” was closely positively related to HPGD in kidney samples of patients with DN (Figure 3F).

Figure 3 The correlation of differentially expressed FA metabolism-related genes with the immune landscape in DN patients. (A) The violin diagram showing the immune cell differences between individuals with or without DN. Red asterisks exhibiting a high proportion of immune cell accumulation in DN, whereas blue asterisks displaying a decrease in DN. (B) The heatmap indicating the correlation among recruited immune cells in DN. (C) The correlation map for the differentially expressed FA metabolism-related genes and recruited immune cells in DN. (D) The association between the expression of PLA2G4A, SOX9, CD36, and C3 with the proportion of T cells gamma delta in DN. (E) The relation between the expression of APOC1 and PTGDS with the proportion of Macrophages M2 in DN. (F) The correlation between the expression of HPGD with the proportion of Neutrophils in DN. *p-value < 0.05; **p-value < 0.01; ***p-value < 0.001; ****p-value < 0.0001. Abbreviation: ns, not significant.

Prediction of Candidate Small Molecular Drugs for Treating FA Metabolism-Based DN Subtypes with Poor Outcomes

Molecular subtyping can be carried out not only to reveal the distinct characteristic of DN but also to facilitate precise therapeutic strategies. As indicated by higher expression of pro-fibrotic and pro-inflammatory factors (Figure 2D), patients in Cluster B exhibited worse outcomes than those in Cluster A. Thus, based on the original therapeutic strategy, supplemented novel drugs for patients in Cluster B were urged to be discovered. First, Up-regulated genes in Cluster B (DN-Cluster B vs Control) were identified according to the cut-off value of adjusted p-value < 0.05 and log2 FC > = 1 (Figure 4A–B). Of note, the overlapped genes between up-regulated genes and genes included in the key module of Cluster B (Additional file: Figure S2G) were imported to the connectivity map (CMAP) database to predict potential small-molecule drugs for further treatment (Figure 4C). By using QUERY CMAP, the top 10 small-molecule drugs, which exhibited the highest negative scores, were identified and considered to be candidate pharmacological therapies for FA metabolism-based Cluster B of patients with DN (Figure 4D). The targeted signaling pathways and chemical structures of these drugs are detailly shown in Figure 4E.

Figure 4 Prediction of the potential small-molecular drugs for treating FA metabolism-based subtypes of DN by CMAP analysis. (A) The volcano plot of the DEGs between control individuals and DN-Cluster B patients. (B) The heatmap of the DEGs between control individuals and DN-Cluster B patients. (C) The Venn diagram identifying 66 key up-regulated genes in DN-Cluster B patients for CMAP analysis. (D) The heatmap showing the top 10 drugs with the highest negative scores in 9 cell lines with CMAP analysis. (E) The chemical structures and detailed description of these top 10 drugs.

Establishment of a Prediction Model Based on FA Metabolism-Related Genes via Machine Learning Algorithms

To screen out the potential diagnostic biomarkers to distinguish patients with DN, the LASSO regression algorithm was employed to reveal 11 potential diagnostic biomarkers out of 39 differentially expressed FA metabolism-related genes (Figure 5A and B). Moreover, using the SVM-RFE algorithm, 7 FA metabolism-related genes were determined as candidate biomarkers (Figure 5C). After intersecting the potential biomarkers extracted from LASSO and SVM-RFE algorithms, only two genes (PRKAR2B and ANXA1) were overlapped via a Venn diagram and were regarded as optimum diagnostic biomarkers (Figure 5D). For better performance in DN prediction, a prediction model was established and visualized as a nomogram based on PRKAR2B and ANXA1 expression by carrying out logistics regression analysis (Figure 5E). Thus, by using the established nomogram, we could convert the renal expression of PRKAR2B and ANXA1 into points, respectively, and the total points were obtained by adding of these two points, corresponding to the risk of DN.

Figure 5 Identification of candidate diagnostic biomarkers and construction of a prediction model based on differentially expressed FA metabolism-related genes by machine learning algorithms. (A and B) LASSO algorithm for identifying potential biomarkers for DN diagnosis. The number of genes (N = 11) was selected as the most suitable parameter for DN prediction. (C) SVM-REF algorithm for screening candidate biomarkers for DN diagnosis. The number of genes (N = 7) was chosen as the optimal argument for DN diagnosis. (D) The Venn diagram showing two overlapped reliable diagnostic biomarkers between LASSO and SVM-REF methods. (E) The prediction model was constructed based on the two diagnostic biomarkers.

Predictive Ability and Performance of the Constructed Nomogram Based on FA Metabolism-Related Genes

The area under the curve (AUC) values of PRKAR2B, ANXA1, and nomogram were 0.95, 0.80, and 0.97, respectively (Figure 6A–C). Moreover, the DN risk prediction nomogram displayed ideal consistency between predictions and actual observations in the DN cohort (Figure 6D). Additionally, the decision curve analysis (DCA) for the nomogram suggested that clinical decision-making based on the nomogram might be beneficial for the prediction of the risk of patients with DN (Figure 6E). Furthermore, the receiver operating characteristic (ROC) curve revealed that the nomogram based on FA metabolism-related genes exhibited a relatively good predictive value of patients with early stages of DN (Figure 6F), implying that the prediction model could also be extensively applied to diagnose the occurrence of early DN.

Figure 6 Verification of the predictive ability and performance of the constructed prediction model for DN risk. (A–C) The ROC curve evaluating the diagnostic efficacy of PRKAR2B (A) and ANXA1 (B) and the constructed prediction model (C) in the risk of DN prediction. (D) The calibration curve of the prediction model for diagnosing DN patients. (E) The DCA for the prediction model. (F) The ROC curve showing the diagnostic performance of the prediction model in the external dataset with the early stage of DN.

The Relevance Between FA Metabolism-Related Genes and Clinical Characteristics of DN

To further elaborate the clinical significance of FA metabolism-related genes in DN, we performed the correlation analysis between the renal expression of PRKAR2B and ANXA1 with clinicopathological features, including serum creatinine (SCr) together with glomerular filtration rate (GFR) using Nephroseq v5 online tool. The results showed that the renal PRKAR2B expression was negatively correlated with the SCr level (Figure 7A) and positively associated with GFR in patients with DN (Figure 7B), suggesting that PRKAR2B may exert renoprotective effects in patients with DN. In addition, high renal ANXA1 expression was associated with the increased SCr level and decreased GFR of DN patients (Figure 7C and D), which revealed that renal ANXA1 might be involved in the progression of DN.

Figure 7 Correlation between renal expression of FA metabolism-related genes and clinical characteristics in DN patients. (A and B) The correlation between renal PRKAR2B expression and the SCr level (A) or GFR (B) in DN. (C and D) The association of renal ANXA1 expression with the SCr level (C) or GFR (D) in DN.

Validation of Data in the Animal Model of DN

The renal expression levels of PRKAR2B and ANXA1 were obtained from the integrated DN microarray datasets, and the results showed that PRKAR2B was significantly decreased, while ANXA1 was increased in patients with DN (Figure 8A). To confirm the results of bioinformatics analysis, we firstly established an experimental STZ/High-fat diet (HFD)-induced diabetic mouse model of DN (Figure 8B). We found that multiple renal functional indicators, including blood glucose, serum creatinine (SCr), blood urea nitrogen (BUN), and urine albumin as well as urine albumin/creatinine (UACR), were significantly elevated, and urine creatinine (UCr) was markedly decreased in the diabetic mice, compared to the control mice (Figure 8C). Moreover, PAS staining showed mesangial matrix expansion and increased mesangial matrix index, and Masson’s trichrome staining displayed excessive collagen deposition in both glomeruli and renal interstitium in diabetic mice, indicating successful establishment of the DN animal model (Figure 8D-E). Next, we assessed the mRNA and protein levels of PRKAR2B and ANXA1 in the DN and control mice. Consistent with the bioinformatics results of DN datasets, the DN animal model displayed similar renal expression patterns of PRKAR2B and ANXA1, in which the mRNA and protein levels of PRKAR2B were significantly decreased, whereas those of ANXA1 were increased in the kidney tissues of the DN mice compared to the control mice (Figures 8F, 9A and B).

Figure 8 Validation of data in the experimental DN animal model. (A) The expression levels of PRKAR2B and ANXA1 in the integrated DN microarray datasets. (B) The flow diagram for establishing an experimental animal model. (C) The blood glucose, serum creatinine, blood urea nitrogen, urine creatinine, and urinary albumin were measured, and the urinary albumin/creatinine (UACR) was calculated. (D) The PAS and Masson’s trichrome staining of the control and DN kidney samples. Bar scale, 50 μm. (E) The quantification of mesangial matrix index and fibrotic area between the control and DN kidney samples. (F) RT-qPCR detected the mRNA levels of PRKAR2B and ANXA1 in the kidney tissues of mice with STZ/HFD-induced DN and in those of the control mice. Data are shown as mean ± SD from groups of six mice. **p-value < 0.01; ***p-value < 0.001; ****p-value < 0.0001.

Figure 9 Detection of the protein levels of PRKAR2B and ANXA1 in the experimental DN animal model. (A) The protein levels of PRKAR2B and ANXA1 in the kidney tissues of mice with STZ/HFD-induced DN and in those of the control mice. (B) Quantitative grey values of protein expression are shown. Data are shown as mean ± SD from groups of six mice. ****p-value < 0.0001.

Discussion

With the development of the second-generation sequencing technology and the standardized establishment of public databases, more and more DN-related transcriptome information has been found and clarified, providing scientific guidance for early diagnosis of diabetes mellitus progression to DN and the exploration of DN pathogenesis. Recent evidence implies that FA metabolism-related genes may participate in the pathogenetic process of DN, but the underlying mechanism of FA metabolism-mediated DN progression has not been fully explained in the existing studies. In this study, we performed an in-depth study on the role of FA metabolism-related genes in DN by integrated bioinformatics analysis and machine learning. A total of 39 differentially expressed FA metabolism-related genes were identified based on analysis of the combined dataset from 54 DN and 57 control kidney samples. Importantly, patients with DN were classified into two subclusters according to the FA metabolism signature. Enrichment analysis of the DEGs and abundance analysis of immune cells between the two subtypes revealed that the significant differences were the immune responses and ECM deposition. Moreover, we also discovered potential drugs for the precise treatment of FA metabolism signature-classified DN subtypes. Furthermore, the FA metabolism-related diagnostic nomogram model was established and evaluated comprehensively, unveiling that this model effectively distinguished patients with DN, regardless of the early or advanced stages of DN. These findings may provide a comprehensive overview of FA metabolism in diabetes-induced kidney injuries and may be beneficial to uncover the underlying mechanisms and precise medication of DN.

Impairment of lipid metabolism, especially FA metabolism, may be involved in kidney diseases (eg, acute/chronic kidney diseases and glomerulonephritis), which has started to attract public attention.^41–44 Inhibiting the oxidative process of FA improves kidney function in DN mice,²⁵ suggesting that dysregulation of FA metabolism may also contribute to the onset and progression of DN. However, the regulatory role of FA metabolism in DN remains obscure. Here, two subtypes of DN were identified based on differentially expressed FA metabolism-related genes. The enrichment analysis based on the DEGs between these two groups suggests that the underlying mechanism of FA metabolism in DN may be closely associated with “inflammatory/responses” and “extracellular matrix organization”. Increasing studies have demonstrated that the recruitment of multiple types of immune cells in the kidney can be observed in different experimental DN animal models^45,46 and clinical biopsy samples at all stages of DN.⁴⁷ Both the innate immune system (eg, macrophages and dendritic cells)⁴⁸ and the adaptive immune system (eg, T cells and B cells) may participate in the progression of DN.⁴⁹ In this study, we noticed that the alteration of the immune landscape was also closely correlated to FA metabolism-related genes-based groupings, as indicated by infiltrated immune cell analysis in the present study. A total of five types of immune cells, including “Dendritic cells resting”, “Plasma cells”, “Treg”, “Macrophages M2”, and “T cells gamma delta”, were differentially recruited in two distinct DN subgroups. Furthermore, the expression of FA metabolism-related genes was associated with the abundance of “Macrophages M2”, “T cells gamma delta”, and “Neutrophils”, demonstrating that FA metabolism was involved in the regulation of immune cell accumulation in patients with DN.

The concept of precise medication has facilitated the classification of patients. Different subtypes may display different pathogenic mechanisms and clinical outcomes. It has been widely accepted that activation of ECM organization leads to the accelerated accumulation of ECM and the progression of renal fibrosis in DN, indicating poor outcomes of kidney function.⁵⁰ This study found the increased expression of pro-fibrotic factors (eg, FN1, COL1A1, and ACTA2) in patients within Cluster B compared to Cluster A. Given the poor clinical outcomes of patients in Cluster B, potential complementary drugs on the basis of typical therapeutic methods for precise treatment need to be explored. During the last few decades, small-molecular drugs targeting a series of diseases have achieved great success.^51,52 Several merits, including high tissue penetration, a tunable half-life, and oral bioavailability, have made small-molecular drugs more effective in conducting precise medication. Using CMAP analysis, several small-molecule drugs for the precise medication of patients in Cluster B were screened based on key genes included in both subsets (up-regulated genes and crucial module genes of Cluster B), whereas further trials are needed for assessing the therapeutic effects of those drugs in ameliorating kidney function of DN patients in Cluster B.

So far, challenges still exist in the diagnosis of DN due to the limited ideal clinical indicators and the heterogeneity of lesions in the kidney.⁵³ Thus, potential biomarkers which exhibit high specificity and sensitivity are still urgently identified. Recently, FA metabolism-related genes have been suggested to be helpful for the risk prediction of various cancers.^54–56 However, the diagnostic performance of FA metabolism-related genes for DN has yet to be explored. In this study, machine learning methods, including LASSO and SVM-REF, were used to explore the candidate FA metabolism-related biomarkers to effectively distinguish DN patients. PRKAR2B and ANXA1 with optimal predictive efficacy were eventually identified. The most noteworthy finding in this study was that a more comprehensive prediction model based on these two biomarkers was established to improve the diagnostic efficacy of DN. Furthermore, an external dataset also validated that the diagnostic model was efficient in predicting the risk of patients with DN of early stages, suggesting the FA metabolism-related prediction model constructed in our study is beneficial for diagnosing patients with DN of different stages. Additionally, significant clinical correlations were found between these two biomarkers and renal function such as SCr level and GFR. Thus, both renal PRKAR2B and ANXA1 have the potential to be biomarkers for predicting the severity of DN.

In summary, our results brought a comprehensive overview of FA metabolism in patients with DN. We demonstrated the potential association between FA metabolism signature and infiltrating immune cells as well as extracellular matrix organization, explored precise medication targeting distinct DN clusters based on FA metabolism-related genes, and identified PRKAR2B/ANXA1 as the candidate biomarkers for patients with DN by LASSO and SVM-REF algorithms. Furthermore, a prediction model combining PRKAR2B and ANXA1 was constructed to diagnose DN of early and advanced stages with high accuracy. Therefore, these findings may elucidate the role of FA metabolism in diabetic kidney injury, providing evidence for better management and therapy of patients with DN.