Comprehensive Analysis Reveals Biomarkers Related to Diabetic Peripheral Neuropathy and Its Molecular Mechanism

Introduction

Globally, there is a concerning trend of the escalating prevalence rates and disease burden associated with diabetes, posing significant challenges to public health systems worldwide.¹ The International Diabetes Federation (IDF) diabetes map, published in 2021, reports that a staggering 536.6 million individuals worldwide afflicted with diabetes, constituting a prevalence of 10.5%.² Diabetic peripheral neuropathy (DPN) is one of the most common complications of diabetes, impacting approximately half of all individuals with diabetes.³ DPN manifests as a chronic, symmetrical, and progressive ailment typified by abnormal sensation or/and pain in its early stages. Patients frequently report lancinating, tingling, and burning sensations. These symptoms are often accompanied by depression, anxiety, and sleep disturbances.⁴ Due to the loss of peripheral nerve function, DPN initially affects the longest nerve fibers, with the condition subsequently progressing from distal to proximal regions.⁵ This progression can result in disability, severely diminishing the patient’s quality of life and imposing substantial physical, emotional, and financial burdens on both the affected individuals and the healthcare system.⁶ Therefore, early and timely diagnosis and treatment of DPN are crucial in improving clinical outcomes and alleviating the economic and emotional burdens on patients. Currently, the diagnosis of DPN relies on nerve conduction studies (NCS), which is considered the gold standard.⁷ However, its time-consuming, labor-intensive, and notably expensive nature makes routine clinical application impractical.⁸ Considering the lack of precise early diagnostic markers for DPN, the existing ambiguity in our understanding of the underlying mechanisms of DPN and the scarcity of effective treatment strategies. There exists a compelling imperative to diligently pursue biomarkers for early detection and a comprehensive examination of the molecular mechanisms underlying these markers is imperative. It will further improve the diagnosis, treatment, and prognosis of this disease.

The Mendelian randomization (MR) methodology holds significant promise in evaluating causal relationships between exposures and outcomes in epidemiology.⁹ By utilizing genetic variants as instrumental variables (IVs) in non-experimental data, MR offers a robust approach to investigating such associations. Since genetic alleles are randomly arranged at conception and the germline genotype cannot be modified by disease, MR analyses are not biased by confounding between exposures and outcomes, reverse causality, or measurement error.¹⁰ Therefore, MR emerges as a compelling approach for meticulously investigating the complex mechanisms involved in the pathogenesis of DPN. Most prior studies have employed MR to identify risk factors for DPN.^11,12 However, fewer investigations focusing on the concurrent utilization of transcriptome analysis and MR in the study of DPN.

This study integrates transcriptome sequencing data of DPN in conjunction with MR analysis to scrutinize the correlation and causality between DEGs and the pathogenesis of DPN. This study will utilize enrichment analysis to investigate the relevant functional pathways and the pathogenesis of these genes. Furthermore, the validation of biomarker expression will be conducted using reverse transcription quantitative polymerase chain reaction (RT-qPCR). These comprehensive integrated methodologies, combining transcriptome sequencing data and MR analysis, are strategically crafted to unravel the intricate molecular mechanisms underlying the pathogenesis of DPN. By meticulously scrutinizing the correlation and causality between DEGs and DPN, these methodologies hold immense promise in unveiling novel therapeutic targets and paving the way for innovative interventions to ameliorate the burden of this debilitating condition.

Material and Methods

Statistical Collection of Data

The GSE95849 and GSE185011 datasets were obtained from the Gene Expression Ontology (GEO) database (https://www.ncbi.nlm.nih.gov/geo/).¹³ The GSE95849 dataset (GPL22448) was selected as the training set, wherein blood samples were obtained from six patients diagnosed with DPN and six control subjects.¹⁴ The GSE185011 dataset (GPL24676) served as the validation set, comprising blood samples obtained from five patients diagnosed with DPN and five normal control subjects.¹⁵ The integrative epidemiology unit open genome-wide association study database (IEU OpenGWAS, https://gwas.mrcieu.ac.uk/)¹⁶ was used to procure expression quantitative trait locus (eQTL) data for differentially expressed DPN-related genes (DE-DPN-RGs) which were employed as exposure factors in the MR analysis. The visit time was October 23, 2023. Simultaneously, a search was conducted by entering the keyword “diabetic neuropathy” in this database to obtain the “finn-b-DM_NEUROPATHY” dataset, utilized as the outcome. The dataset comprised 1,415 disease samples and 162,201 control samples, all sourced from European populations. The number of single nucleotide polymorphisms (SNPs) was 16,380,195.

Analysis of Differentially Expressed Genes (DEGs)

Gene expression data were extracted from DPN and normal control samples in the GSE95849 dataset. DEGs1 were then identified through screening using the “limma” (version 3.46.0)¹⁷ with |log₂fold-change (FC)| > 1, adj. p < 0.05.¹⁸ Building upon this foundation, DEGs1 were visualized via the “ggplot2” (version 3.4.4)¹⁹ and the “pheatmap” (version 1.0.12) packages.²⁰ In addition, DEGs2 between DPNs and controls in the GSE185011 dataset were identified and visualized using the same methodology. Following this, the sets of DEGs identified in the GSE95849 and GSE185011 datasets were intersected to obtain the subset of differentially expressed DPN-related genes (DE-DPN-RGs).

Functional Enrichment Analysis of DE-DPN-RGs and Construction of a Protein-Protein Interaction (PPI) Network

The subset of differentially expressed DE-DPN-RGs underwent Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis using the “clusterProfiler” package (version 4.4.4), with a significance threshold of P < 0.05.²¹ This analysis aimed to elucidate the biological functions and signaling pathways associated with the DE-DPN-RGs. Additionally, interactions among the DE-DPN-RGs at the protein level were investigated. The PPI network of the DE-DPN-RGs was constructed using the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://string-db.org),²² with an interaction score threshold greater than 0.7, and visualized using “Cytoscape” software (version 3.8.2).²³

Selection of Instrumental Variables (IVs)

The following three assumptions were followed throughout the MR analysis: 1. The independence assumption (the causality between DE-DPN-RGs and DPN is independent of confounding factors); 2. The association assumption (there exists a robust correlation between IVs and DE-DPN-RGs); 3. The exclusivity assumption (DE-DPN-RGs represent the sole pathway through which genetic variation influences DPN).²⁴ Subsequently, IVs satisfying all three aforementioned assumptions were identified through screening. Initially, the extract_instruments function of the “TwoSampleMR” package (version 0.5.7).²⁵ It was employed to retrieve the differentially expressed DE-DPN-RGs and filter the IVs significantly correlated with the DE-DPN-RGs (p < 5×10⁻⁸). Subsequently, single nucleotide polymorphisms (SNPs) in linkage disequilibrium (LD) were excluded by the Clump function (clump = TRUE, r² = 0.001, kb = 10,000). After that, the number of SNPs was greater than 3 and the SNPs with F-values greater than 10 were reserved. The calculation formula for F values was as follows. In the formula, R² represents the cumulative explanatory variance of the selected SNPs, N represents the number of the samples.

Following this, the extract_outcome_data function was utilized to retrieve the outcome data, and the IVs that were not relevant to the outcome was retained. SNPs meeting stringent screening criteria were then employed for subsequent MR analysis. Moreover, the harmonise_data function was applied to standardize effect alleles and effect sizes, ensuring consistency in allele representation between exposure factors and outcomes. Similarly, SNPs for multivariable Mendelian randomization (MVMR) analysis were selected using the aforementioned methodology.

Univariable MR (UVMR) Analysis

The univariable UVMR analysis of DE-DPN-RGs (exposure factors) and DPN (outcome) were performed to investigate the causal relationship between them. In UVMR analysis, inverse variance weighted (IVW) method²⁶ was employed as the primary approach, supplemented with the weighted median (WM),²⁷ MR Egger,²⁸ simple mode,²⁹ and weighted mode³⁰ methods. Notably, a p-value less than 0.05 from the IVW method indicated a significantly causal relationship between DE-DPN-RGs and DPN. Additionally, odds ratios (ORs) greater than 1 indicated that DE-DPN-RGs were risk factors for DPN, while ORs less than 1 indicated that DE-DPN-RGs were protective factors for DPN. Furthermore, the UVMR results were visualized using scatter plots, forest plots, and funnel plots.

Sensitivity Analysis

Subsequently, the reliability of the UVMR results was verified through sensitivity analyses. Firstly, a heterogeneity test was used to assess the presence of heterogeneity between the two datasets via the “TwoSampleMR”. If the p-value is greater than 0.05, there was no heterogeneity. If the results had heterogeneity (P < 0.05), random effects IVW was adopted in the MR, otherwise fixed effects IVW would be adopted.Secondly, the mr_pleiotropy_test function was utilized to assess whether there was horizontal pleiotropy between DE-DPN-RGs and the DPN dataset. In addition, the MR-PRESSO was used to correct for horizontal pleiotropy between them by removing outliers. If the p-value is greater than 0.05, it meant that there are no confounders and the results were trustworthy. The SNPs were then eliminated one by one using the leave-one-out (LOO) test to determine the effect of the remaining SNPs on the results.

MVMR Analysis

MVMR was further performed on the exposure factors screened by the UVMR analysis to determine if the exposure factors had a direct or indirect effect on DPN. Allelic effects and effect sizes were reconciled for DE-DPN-RGs and DPN data using the mv_harmonise_data function of “TwoSampleMR”. Unnecessary exposure factors were then removed using the mv_lasso_feature_selection function. Finally, the Steiger test was done to rule out the effect of reverse causation.

Identification of Biomarkers for DPN

Genes with a significantly causal relationship with DPN were screened as candidate biomarkers for DPN from DE-DPN-RGs based on UVMR analysis. In addition, the expression of candidate biomarkers in the GSE95849 and GSE185011 datasets was examined and validated. Ultimately, candidate biomarkers with significant expression differences between DPN and controls and consistent expression trends in both datasets could be used as biomarkers for DPN.

Functional Analysis of Biomarkers

In the GSE95849 dataset, the samples were divided into high and low expression groups based on the expression of biomarkers, and the genes between them were analyzed for differences, and |log₂FC| of each gene was calculated and ranked from largest to smallest. The biomarkers were then subjected to gene set enrichment analysis (GSEA)³¹ using the “clusterProfiler” (version 4.7.1.003).³² The reference gene set was the background gene set “c2.all.v7.2.symbols.gmt” from the Molecular Signatures Database (MSigDB)³³ (P < 0.05). To explore differential biological pathways between high and low expression groups of biomarkers, gene set variation analysis (GSVA) was performed on different expression groups based on the “c2.all.v7.2.symbols.gmt” (P < 0.05). In addition, the co-expression network of biomarkers and other functionally similar genes was explored using GeneMANIA (http://genemania.org/). In addition, functional similarities between biomarkers were analyzed with the help of the “GOSemSim” (version 2.24.0)³⁴ (functional similarity score > 0.5).

Subcellular Localization and N6-Methyladenosine Binding Site Prediction

The localization of the biomarkers in subcells was predicted with the help of the mRNALocater online website (http://bio-bigdata.cn/mRNALocater).³⁵ The localization scores of the biomarkers in different subcells were calculated. Next, the occupancy of the biomarkers in different subcells was plotted based on this score. To explore whether N6-methyladenosine (m⁶A) affected the translational stability of biomarkers, the biomarker m⁶A binding sites were predicted. First, genome sequence for the biomarkers were obtained at National Center for Biotechnology Information (NCBI, https://www.ncbi.nlm.nih.gov/).³⁶ Then the m⁶A binding sites of the biomarkers were predicted from the sequence-based RNA adenosine methylation site predictor database (SRAMP, http://www.cuilab.cn/sramp).³⁷

Construction of Biomarker Regulatory Networks

In order to investigate the regulatory mechanisms of biomarkers at the molecular level, the miRDB database (http://mirdb.org)³⁸ and TargetScan database (http://www.targetscan.org)³⁹ were used to predict the upstream miRNAs corresponding to biomarkers. And miRNAs were then obtained by taking the intersection of the predictions from the two databases. After that, the upstream lncRNAs of the miRNAs were predicted using the Starbase database (https://starbase.sysu.edu.cn/).⁴⁰ In addition, a competing endogenous RNA (ceRNA) regulatory network was constructed using the “Cytoscape”. Potential transcription factors (TFs) for biomarkers were predicted in the Human Tranion Factor Targets database (hTFtarget, http://bioinfo.life.hust.edu.cn/hTFtarget/).⁴¹ Next, the TF-mRNA-miRNA network was mapped using the “Cytoscape”.

Prediction of Drugs for DPN-Targeted Therapies

The Drug Signatures Database (DsigDB, http://tanlab.ucdenver.edu/),⁴² an online database, was used to predict the targeted drugs associated with the biomarkers, and the targeting relationships between the biomarkers and drugs were demonstrated using “Cytoscape”.

Expression Analysis of Biomarkers

For validating the expression of biomarkers, reverse transcription quantitative polymerase chain reaction (RT-qPCR) was performed on 5 DPN and 5 control samples from peripheral blood mononuclear cell (PBMC). The blood samples used in this study were obtained from the Second Affiliated Hospital of Zunyi Medical University, with sample collection occurring between December 2023 and April 2024. All subjects provided informed consent prior to participation. The collection process adhered to the ethical standards set by the Medical Ethics Committee of Zunyi Medical University. The amplification conditions for RT-qPCR were 40 cycles with 1 minute at 95°C, 20 seconds at 95°C, 20 seconds at 55°C, and 30 seconds at 72°C. The qPCR primers were listed in Table S1 with GAPDH as reference gene. The relative expression levels of biomarkers were calculated using the 2^−ΔΔCT method.

Statistical Analysis

The R programming language (version 4.2.1) was engaged in statistical analysis. Wilcoxon test was utilized to examine the discrepancies between 2 groups. P < 0.05 was considered notably significant. If not specified, adj. P<0.05 or P<0.05 indicates significant differences.

Results

The Biological Functions of 124 DE-DPN-RGs

Among 5,340 DEGs identified between DPN and controls in the GSE95849 dataset, 3,746 DEGs exhibited upregulation, while 1,594 DEGs demonstrated downregulation. The distribution of these DEGs1 was visualised using a volcano and a heatmap (Figure 1A and B). A total of 896 DEGs2 were identified between DPN and controls in the GSE185011 dataset, with 823 DEGs2 were up-regulated and 73 DEGs2 were down-regulated (Figure 1C and D). The intersection of the 3,746 DEGs from the GSE95849 dataset with the 896 DEGs from the GSE185011 dataset yielded a total of 124 DE-DPN-RGs (Figure 1E). Additionally, DE-DPN-RGs were significantly enriched in several biological functions, including the “intrinsic apoptotic signaling pathway”, “cell chemotaxis”, “regulation of neuron death”, “positive regulation of cytokine production”, “positive regulation of interleukin-6 production”, and “Toll-like receptor binding” (P < 0.05) (Figure 1F). In the KEGG analysis, DE-DPN-RGs were significantly enriched in the “TNF signaling pathway”, “Toll-like receptor signaling pathway”, “Antifolate resistance”, and “African trypanosomiasis” (P < 0.005) (Figure 1G). Furthermore, a robust linkage among these 124 DE-DPN-RGs was demonstrated in the PPI network, with prominent interactions observed among genes TNF, HIF1A, NOD2, and SAMD9L (Figure 1H).

Figure 1 Results of differential expression analysis. (A) The volcano plot of differential genes in the GSE95849 dataset, with blue representing down-regulated genes and purple representing up-regulated genes; (B) The heatmap plot of differential genes in the GSE95849 dataset; (C)The volcano plot of differential gene expression in the GSE185011 dataset; (D) The heatmap plot of differential genes in the GSE185011 dataset; (E) The Venn diagrams for DEG1 and DEG2; (F) Results of GO enrichment analysis (P < 0.05); (G) Results of KEGG enrichment analysis (P < 0.005); (H) Results of PPI network.

TNF, OSBPL8, IER3, SLC16A3, CREB5, and LRWD1 Were Significantly and Causally Associated with DPN

A total of six DE-DPN-RGs causally associated with DPN were screened as candidate biomarkers for DPN. Among them, TNF was a protective factor for DPN (OR < 1, P < 0.05). OSBPL8, IER3, SLC16A3, CREB5, and LRWD1 were risk factors for DPN (OR > 1, P < 0.05) (Figure 2A). The scatter plot further demonstrated this result. TNF had a negative slope and was a protective factor for DPN, whereas OSBPL8, IER3, SLC16A3, CREB5, and LRWD1 had positive slopes and were risk factors for DPN (Figure S1). The results based on the IVW algorithm showed that the distribution of SNP loci for the diagnostic efficacy of DPN in the data of TNF was all less than zero in the forest plot, which further suggested that TNF was protective factor for DPN. In contrast, the distribution of SNP loci diagnostic of DPN efficacy in the OSBPL8, IER3, SLC16A3, CREB5, and LRWD1 data were all greater than zero, which further suggested that all of these genes were risk factors for DPN (Figure S2). SNPs were uniformly and randomly distributed on both sides of the IVW line, suggesting that the analyses followed Mendel’s second law (Figure S3).

Figure 2 The forest map of candidate biomarkers. (A) Six DE-DPN-RGs exhibit causal associations with DPN; (B) Results of multifactorial Mendelian randomisation analysis.

The results of the heterogeneity test indicated that there was no heterogeneity between the DE-DPN-RGs (OSBPL8, SLC16A3, CREB5, LRWD1, and TNF) dataset and the DPN dataset (p > 0.05). Although the p-value of IER3 was less than zero, the heterogeneity did not affect the results too much due to the IVW algorithm used (Table 1). Moreover, the results of the horizontal polytropy test proved that there was no influence of confounding factors in the analysis of this study (p > 0.05) (Table 2). Also, MR-PRESSO results confirmed the absence of horizontal pleiotropy between them (Table 3). Then, in the causal relationship between these six candidate biomarkers and DPN did not change significantly by excluding a SNP, which confirmed the stability of the result (Figure S4). Additionally, the six candidate biomarkers passed the Steiger test (p < 0.0001), confirming the reliability of the forward analysis of causality (Table 4). Finally, the results of the MVMR analysis showed that TNF remained a protective factor for DPN, and OSBPL8, IER3, SLC16A3, CREB5, and LRWD1 remained risk factors for DPN (Figure 2B).

Table 1 Results of MR Heterogeneity Test Analysis

Table 2 Results of MR Level Multiple Validity Test Analysis

Table 3 Results of the Analysis of the MR-PRESSO Horizontal Multiple Validity Test

Table 4 Reverse Causality Analysis

CREB5, LRWD1, OSBPL8, and SLC16A3 Were Biomarkers of DPN

Validation in the GSE185011 dataset showed that the expression levels of CREB5, LRWD1, OSBPL8, and SLC16A3 were significantly different between DPN and controls, with consistent trends in the GSE95849 dataset, and that they were upregulated in DPN (Figure 3A and B). Therefore, these four genes were credited as biomarkers of DPN. Subsequently, the GSEA results demonstrated the signaling pathways in which the biomarkers were involved. GSEA results showed that the signaling pathways co-enriched by these four biomarkers were “metabolic reaction groups of amino acids and derivatives”, “rRNA processing”, and “translation” (Figure 3C–F).

Figure 3 The results of expression validation analyses and GSEA enrichment analyses. The expression results of biomarkers in GSE185011 dataset (A) andGSE95849 (B); *: P < 0.05, **: P < 0.01, ***: P < 0.001, ns, no significance. The results of GSEA analysis of four biomarkers (C–F).

Functions and Functional Associations of Biomarkers

GSVA results showed significantly different biological pathways. Among them, “reactive protein NF-κb activation through fadd rip 1 pathway bmediated by caspase 8 and 10”, “biocarta sppa pathway” and “reactome clec7a dectin 1 induces nfat activation” were significantly highly expressed in the high expression group. For example, signaling pathways such as “reactome phospholipase c mediated cascade FGFR4”, and “NF-κb is activated and signals survival”, “PI3K cascade FGFR3”, were significantly highly expressed in the low expression group (Figure 4A–D). The GeneMANIA results showed that there were various relationships between biomarkers and other genes such as co-expression, the same location, physical interactions, genetic interactions, sharing protein structural domains, or participating in the same pathway. Among them, LRWD1 interacted most closely with ORC5 and ORC2 (Figure 4E). Furthermore, the results of the Friends analysis demonstrated a low functional similarity between the four biomarkers (functional similarity score < 0.5) (Figure 4F).

Figure 4 The GSVA enrichment analysis and GeneMANIA database analysis. (A–D) The Heatmap of GSVA analysis of four biomarkers; (E) The GeneMANIA analysis results; (F) The Results of Fhorends analysis.

Subcellular Localization of Biomarkers and Where Biomarkers Undergo m6A Methylation

SLC16A3, CREB5 and LRWD1 were detected in the cytoplasmic region, whereas OSBPL8 was detected as up-regulated in the extracellular region. SLC16A3 was also detected in the endoplasmic reticulum (Figure 5A). Meanwhile, the predicted locations of biomarkers susceptible to m⁶A methylation showed that CREB5 was more susceptible to m⁶A methylation at positions 0 to 1000 of the base (Figure 5B). LRWD1 was prone to m⁶A methylation at about 1000,1500, and 2000 base positions (Figure 5C). SLC16A3 was found to be susceptible to m⁶A methylation at around 1500 to 2000 base positions (Figure 5D). OSBPL8 was susceptible to m⁶A methylation at around 60 base positions (Figure 5E).

Figure 5 The results of subcellular localisation and m6A binding site prediction of key genes. (A) Subcellular localisation results for biomarkers; The results of methylation analysis of four biomarkers, (B) CREB5 (C) LRWD1 (D) SLC16A3 (E) OSBPL8.

Molecular Regulatory Network of DPN and Its Potential Target Drugs

A total of 54 miRNAs were obtained by taking the intersection of the two database predictions. A total of 759 lncRNAs were predicted based on 54 miRNAs. Based on this, a ceRNA network was constructed with OSBPL8 and CREB5 (Figure 6A). Meanwhile, CREB5 and OSBPL8 predicted 28 TFs and 20 TFs, respectively, and a TF-mRNA-miRNA regulatory network was constructed by combining the predicted miRNAs (Figure 6B). Where PBX3, ETS1, EP300, CREBBP, and BRD4 were TFs predicted by CREB5 and OSBPL8. In addition, there are 43 targeted drugs available for DPN treatment (Figure 6C). Examples included rosuvastatin, nitrofurfurylhydrazine, and neostigmine bromide.

Figure 6 The ceRNA regulatory network analysis, TF-mRNA-miRNA network construction and drug prediction results. (A) A Map of a ceRNA network constructed by OSBPL8 and CREB5; (B) TF-mRNA-miRNA regulatory network; (C) A targeted drug prediction map for DPN treatment based on OSBPL8 and CREB5.

All Four Biomarkers Were Highly Expressed in the DPN Group by RT-qPCR

As verified by RT-qPCR, the expression levels of CREB5, LRWD1, OSBPL8, and SLC16A3 were higher in the DPN group, which was consistent with the trend in the GSE185011 and GSE95849 datasets, with only OSBPL8 showing a non-significant difference between groups (Figure 7).

Figure 7 The expression results verified by RT-qPCR. *: P < 0.05, **: P < 0.01, ns, no significance.

Discussion

As one of the most prevalent complications of diabetes, DPN significantly impairs patients’ quality of life.⁴³ The therapeutic efficacy of DPN is often suboptimal, and the treatment duration is typically prolonged. Therefore, early diagnosis plays a pivotal role in facilitating timely treatment initiation, which in turn contributes to achieving favorable treatment outcomes, and alleviating the associated economic burden. However, the absence of precise indicators for early DPN detection has hindered the development of accurate diagnostic modalities.⁸ In this study, we identified six candidate biomarkers, TNF, OSBPL8, IER3, SLC16A3, CREB5, LRWD1 through transcriptomics and MR techniques. We found that TNF is a protective factor for DPN. Conversely, the remaining five genes were identified as risk factors for the development of DPN. To enhance the stability and reliability of our findings, we conducted expression validation analyses for these six candidate biomarkers. Among them, the four candidate biomarkers with significant expression and the same trend between the training and validation sets were recorded as the biomarkers: OSBPL8, SLC16A3, CREB5, and LRWD1. Additionally, elucidating the molecular mechanisms holds paramount importance for advancing clinical research on DPN and enhancing our understanding of its pathogenesis.

OSBPL8, belonging to the oxysterol-binding protein-like (OSBPL) family, serves as a crucial intracellular lipid binding and transport protein involved in lipid transportation and cholesterol level regulation.⁴⁴ OSBP-related protein 8 (ORP8), encoded by OSBPL8, is targeted for silencing by miR-143, a micro-RNA associated with obesity, diabetes and cancer.⁴⁵ Decreased expression of ORP8 in cultured hepatocytes has been shown to impairaffects insulin-induced AKT activation. Furthermore, induced transgenic overexpression of miR-143 impairs insulin-stimulated AKT activation and glucose homeostasis, suggesting that the miR-143-ORP8 pathway may represent a potential target for obesity-associated diabetes mellitus.⁴⁶ Specifically, this study demonstrates a causal relationship between OSBPL8 and DPN. We speculate that OSBPL8 may affect DPN through the miR-143-ORP8 pathway.

SLC16A3 (MCT4) belongs to the monocarboxylate transporter (MCT) family, facilitating the transport of pH-dependent monocarboxylates such as lactic acid and pyruvic acid.⁴⁷ Recent investigations indicate that lactic acid enhances neuronal mitochondrial metabolism and energy generation in peripheral nerves.⁴⁸ DPN apparently involves metabolic dysfunction and energy depletion of multiple cells within the peripheral nervous system, implying a robust association between SLC16A3 and DPN. Moreover, the spatial distribution of MCT within the peripheral nervous system (PNS) hints at its potential involvement in lactic acid transport and myelination.⁴⁹ Notably, a study revealed that the experimental exacerbation of DPN was associated with a reduction in monocarboxylate transporter-1 (MCT1).⁵⁰

CREB5 is a member of the cAMP-response element binding (CREB) protein gene family and is classified under the ATF2 subfamily. This subgroup demonstrates a heightened affinity for cAMP-response elements (CREs) and participates in diverse cellular functions such as survival, proliferation, and glucose metabolism.⁵¹ Recent studies have suggested a potential protective role for activated CREB5 in shielding diabetic proximal renal tubule cells from apoptosis.⁵² Additionally, CREB5 has been associated with the AMPK signaling pathway,⁵³ which is speculated to be pivotal in the pathogenesis of DPN. AMPK-mediated mitochondrial impairment in Schwann cells may instigate demyelination and axonal atrophy, both implicated in the development of the disease.⁵⁴

LRWD1 is an origin recognition complex (ORC)-associated protein that plays a key role in heterochromatin replication in mammalian cells by recruiting promoter ORCs.⁵⁵ Few previous studies have associated it with diabetes-related diseases, and this study is the first to find an association between LRWD1 and DPN.

Based on our differential analysis, four candidate biomarkers, CREB5, LRWD1, OSBPL8, and SLC16A3, exhibited significant differential upregulation. Moreover, beyond the the causal relationships found by MR, the consistency of their expression trends was verified in both the training and validation sets. Given these insights, these genes emerge as promising targets for therapeutic intervention in DPN. Further investigation into their molecular mechanisms is imperative to enhance our understanding and unlock their therapeutic potential.

Of particular interest in our GO analysis was toll-like receptor binding. Toll-like receptors (TLRs) are pivotal proteins in the innate immune response, initiating a signaling cascade upon interaction with microbial pathogen-expressed TLR ligands. This cascade leads to the production of cytokine and the initiation of adaptive immune responses. Dysregulated immune mechanisms play a crucial role in the pathogenesis of DPN.⁵⁶ Among the TLR, TLR4 plays a crucial role in inflammatory diseases.⁵⁷ DPN is increasingly recognized as a chronic inflammatory condition wherein TLR4 expression in dorsal root ganglia can be induced by high glucose, free fatty acids, and angiotensin II. This induction activates downstream signaling pathways of inflammation, leading to the activation of NF-κB and the subsequent production of inflammatory factors such as TNF-α, prostaglandins, and activated cyclooxygenase (COX), thereby instigating neuroinflammation in DPN.⁵⁸ Notably, the Toll-like receptor signaling pathway was similarly enriched in our KEGG results. Moreover, TLR4 has been identified as a potentially sensitive diagnostic biomarker for DPN in patients with type 2 diabetes.⁵⁷ This evidence underscores the significant role of Toll-like receptors in the progression of DPN and suggests that DEGs may influence DPN through the Toll-like receptor signaling pathway.

Our biomarkers underwent further analysis using GSEA to elucidate their molecular mechanisms. They were notably enriched in pathways related to metabolic reactions of amino acids and derivatives, rRNA processing, and translation. While the metabolism of amino acids and their derivatives in DPN has been scarcely described, studies have highlighted their significance in relevant pathways in diabetic foot ulcers (DFU). Furthermore, literature has underscored a significant association between DPN and DFU,⁵⁹ with DPN serving as a significant initiating risk factor for DFU development.⁶⁰ Although a causal relationship between the two variables was not established in our study, it’s worth noting their interconnectedness. Moreover, GSVA analyses revealed enrichment of the NF-κB pathway, consistent with the potential for neuroinflammation in DPN via the production of inflammatory factors. This further underscores the importance of the signaling pathways involving key genes in DPN progression. Additionally, CREB5 is markedly activated in proximal renal tubular epithelial cells of diabetic patients⁵² and has been found to regulate the MET signaling pathway crucial for colorectal cancer metastasis.⁶¹ OSBPL8 overexpression has been observed in cholangiocarcinomas.⁶² SLC16A3 is abnormally upregulated at the plasma membrane of cardiomyocytes from patients with type 2 diabetes⁶³ and its role in intracellular pH regulation and lactate-based metabolism has been explored in various cancers, including colorectal, breast, and lung cancers.⁶⁴ In summary, these biomarkers are not only implicated in the pathogenesis of DPN and diabetes-related diseases but also play roles in cancer, providing potential avenues for the treatment of related cancers.

In order to delve into the genes linked with biomarkers functions and elucidate the functional correlations among different genes, we conducted GeneMANIA analysis alongside functional correlation analysis. GeneMANIA analysis revealed that LRWD1 interacts most closely with ORC5 and ORC2. LRWD1 interacts directly with the N-terminus of ORC2, and the stability of LRWD1 depends on its binding to ORC2.⁶⁵ The interaction of LRWD1 with ORC5 has been rarely mentioned in previous studies. Notably, the functional similarity among these four biomarkers was observed to be modest. Furthermore, we undertook subcellular localization analysis and predicted m6A binding sites for the biomarkers to gain deeper insights into their roles. The subcellular localization provides important information about the function and activity of biomarkers in the cell. DNA methylation of the SLC16A3 promoter regulates MCT4 expression and affects clinical outcomes in renal cancer.⁶⁶ Hypomethylation of CREB5 upregulates CREB5, thereby altering trophoblast cell function.⁶⁷ Thus, m6A binding site prediction results indicate that m6A methylation of these genes at different base positions may have an impact on their function and regulation. It provides important guidance and insights for further studies in the future.

To delve deeper into the functions and regulatory mechanisms of the biomarkers, we analyzed competing endogenous RNA (ceRNA) regulatory networks. In this regard, we identified 54 miRNA target genes and predicted 759 lncRNA target genes based on these miRNAs to construct a lncRNA-miRNA-key gene regulatory network. Notably, HSA-miR-23a-3p, a member of a major miRNA family involved in regulating T-cell immune responses, exhibited significant reduction in diabetic patients. Inhibition of the NF-κB pathway by HSA-miR-23a-3p has been shown to suppress the expression of pro-inflammatory molecules and angiogenic potential, potentially ameliorating neuroinflammation in DPN.⁶⁸ Additionally, miR-23a has been implicated in inhibiting high glucose-induced inflammation in microglia by down-regulating PDE4B, attenuating DPN.⁶⁹ Moreover, in individuals with type 1 diabetes and long-term metabolic syndrome, HSA-miR-495-3p is downregulated, and previous studies have linked this miRNA to immune and inflammatory responses.⁷⁰ Notably, in patients with DFU, tissue-specific expression of miR-23a and miR-23b was low, while miR-23c was elevated.⁷¹ Considering the association of DPN with DFU, we posit that HSA-miR-23a-3p, HSA-miR-495-3p, HSA-miR-23b-3p and HSA-miR-23c could potentially be ceRNA pathways implicated in DPN regulation. However, further experiments are warranted to confirm this hypothesis.

Additionally, we constructed a transcription factor (TF) regulatory network of key genes to explore the relationships between transcription factors, miRNAs, and key genes. Predicted TFs by OSBPL8 and CREB5 included PBX3, ETS1, EP300, CREBBP, and BRD4. Among them, CREBBP and EP300 have been found to be significant for DPN diagnosis, potentially affecting DPN by influencing lipid metabolism toxicity or inflammatory processes triggered by inflammatory stimuli.⁷² Furthermore, BRD4 plays a pivotal role in the pathogenesis of high glucose-induced cardiomyocyte hypertrophy through the AKT pathway.⁷³ Based on existing studies, it’s evident that some predicted transcription factors are indeed involved in DPN regulation, shedding light on the pathogenesis of DPN. This offers new perspectives and clues for a comprehensive understanding of DPN pathogenesis.

In addition, we conducted drug predictions for the four biomarkers, yielding a total of 43 potential drugs. Among these, examples include Rosuvastatin and Neostigmine Bromide. Studies have recognized that Rosuvastatin may alleviate DPN through its anti-fibrotic synergy.⁷⁴ Neostigmine bromide functions by inhibiting the enzymatic cleavage of acetylcholine, thereby reversibly inhibiting acetylcholinesterase. Consequently, it can persist in cholinergic nerve endings for an extended period, eliciting an excitatory effect on smooth and skeletal muscles. Furthermore, we substantiated biomarker expression through RT-qPCR, with outcomes aligning consistently with the observed trends across the two datasets employed in our study. Four biomarkers were highly expressed in the DPN group, with only OSBPL8 showing a non-significant difference between groups. This observation may be attributed to the constrained sample size utilized for RT-qPCR validation.

It is noteworthy that the present study has certain limitations that warrant attention in subsequent research. Our analysis was conducted using retrospective data derived from public datasets, which may introduce case selection bias and limit the generalizability of our findings. Consequently, the diagnostic and predictive utility of biomarkers in DPN requires further validation through large-scale clinical trials. Additionally, the regulatory mechanisms of these biomarkers in DPN need to be further elucidated through functional and mechanistic studies.

Conclusion

In conclusion, our study leveraged transcriptomic analysis to pinpoint genes exhibiting significant differential expression between disease and control groups. Subsequently, we employed both univariate and multivariate MR studies. This allowed us to elucidate the causal relationships between DE-DPN-RGs and DPN, yielding more reliable findings and and the identification of four potential biomarkers related to DPN: SLC16A3, CREB5, OSBPL8 and LRWD1. Moving forward, we intend to closely monitor the progress of these biomarkers and targeted drugs. Our study possesses several noteworthy limitations that warrant consideration. First, the clinical utility of biomarkers necessitates validation across diverse patient populations. Moreover, the targeted drugs identified in our study require thorough validation through extensive clinical trials before obtaining approval for clinical application. Ultimately, our findings offer new insights into the prevention and treatment of DPN, paving the way for novel strategies and directions in early diagnosis.