Identification of novel therapeutic targets for diabetic neuropathy through integrated proteomics and transcriptomics approaches

Abstract

BACKGROUND

Diabetic neuropathy (DN) is a progressive disorder with limited effective treatment options.

AIM

To identify potential therapeutic targets for DN by integrating plasma proteomic and transcriptomic data.

METHODS

A comprehensive analytical framework was developed to identify multi-omics biomarkers of DN. Protein-protein interaction network and Gene Ontology analyses were performed to explore the biological functions of biomarkers. Tier 1 target proteins were further analyzed. Candidate drug prediction and molecular docking studies were conducted to identify potential treatments while assessing the side effects of key target proteins. The mediation of immune cells in the association between proteins and DN was examined through two-step network Mendelian randomization (MR) analysis.

RESULTS

Nine DN-associated proteins were identified by analyzing protein quantitative trait loci from extensive genome-wide association study data. BTN3A1 and MICB were confirmed using MR, summary data-based MR, and colocalization analyses. Of the nine, HSPA1B, PSMB9, BTN3A1, SCGN, NOTUM, and MICB showed negative associations with DN, whereas WARS, BRD2, and CSNK2B were positive. Gene Ontology analysis indicated enrichment in inflammatory response and neuronal injury pathways. BTN3A1 and MICB were identified as Tier 1 targets. Drug prediction and molecular docking analyses indicated cyclosporin A as a potential therapeutic candidate. Two-step network MR analysis showed that MICB mediated DN through human leukocyte antigen-DR++ monocytes. These integrated findings point to an immune-mediated mechanism with translational potential and nominate BTN3A1 and MICB for focused functional validation.

CONCLUSION

Our integrated multi-omics approach identified two promising therapeutic targets for DN, laying the groundwork for new treatment strategies and enhancing our understanding of MICB’s role in DN.

Key Words: Diabetic neuropathy; Proteome; Transcriptome; Therapeutic targets; Mendelian randomization; Multi-omics

Core Tip: This integrated multi-omics study identified nine diabetic neuropathy (DN)-associated plasma proteins, validated across independent datasets. BTN3A1 and MICB emerged as top-tier therapeutic targets, robustly confirmed via Mendelian randomization (MR), summary-data-based MR, and colocalization analyses. Functional enrichment links these proteins to inflammation and neuronal injury. Crucially, cyclosporine A was predicted as a potential repurposed drug candidate targeting MICB. Furthermore, two-step network MR revealed MICB mediates DN risk specifically through human leukocyte antigen-DR++ monocyte abundance, providing a novel mechanistic immune pathway. This work unveils BTN3A1 and MICB as promising targets and elucidates MICB’s immune-mediated role in DN pathogenesis.

INTRODUCTION

Diabetic neuropathy (DN) encompasses a spectrum of neuropathic conditions characterized by various clinical presentations[1]. DN is highly prevalent and can significantly affect patients by elevating their risk of falls, inducing pain, and diminishing their quality of life (QOL). By 2045, it is estimated that 783 million adults worldwide will be affected by diabetes[2]. Despite the severity and widespread occurrence of this complication, effective clinical treatment remains elusive. Hence, innovative strategies are imperative to enhance patient survival and improve QOL beyond symptom management. Identifying and understanding the biomarkers involved in DN pathways are crucial for identifying potential therapeutic targets.

Large-scale genome-wide association studies (GWAS) have paved the way for omics research[3-5]. Mendelian randomization (MR) can be applied to infer causal relationships between exposure and outcome[6,7]. This approach facilitates the validation of hypotheses that lack strong theoretical support. Previous studies and reviews have suggested a potential correlation between vitamin D levels and DN. Nonetheless, definitive evidence remains lacking because of small sample sizes and the inherent challenges of conducting large-scale clinical investigations. The integration of GWAS and MR approaches provides an opportunity to leverage genetic data to robustly assess and potentially confirm the causal relationship between vitamin D and DN[8]. Although systemic inflammation has long been recognized as a pivotal mechanism in DN, the specific inflammatory factors involved have not been well elucidated. A recent MR analysis identified several inflammatory factors as significant contributors to DN development[9]. However, establishing a causal relationship is insufficient, and deeper mechanistic insights are essential for the identification and validation of potential therapeutic targets for DN.

Advancements in oligonucleotide and immunoassay platforms, such as SomaScan and Olink, have enabled the detection of various proteins, widely ranging from 1000 to 7000. Large-scale plasma proteomic GWAS datasets involving substantial sample sizes, including studies on 35559 individuals from Iceland and 54306 participants from the United Kingdom Biobank (UKB), have been documented[10-12]. These datasets offer an opportunity to explore potential drug targets for DN treatment. Previous studies have employed plasma proteomics to investigate drug targets for conditions such as DN[13], chronic kidney disease[14], diabetes[15], Alzheimer’s disease[16], migraines[17], and cancer[18]; nonetheless, there is a dearth of research focusing on proteomic-based drug targets for DN. Therefore, our objective was to identify key therapeutic targets for DN through large-scale plasma proteomic analysis and to further validate the corresponding coding genes of these targets using transcriptome-based GWAS datasets such as eQTLGen[19] and Genotype-Tissue Expression (GTEx)[20].

The fusion of genomics, transcriptomics, and proteomics holds significant promise for determining whether selected biomarkers are involved in causal pathways. In this study, we integrated extensive GWAS datasets pertaining to DN, blood proteomics, and tissue-specific gene expression. Using a comprehensive analytical approach, we extensively explored fundamental biological processes associated with DN. Furthermore, for translational applications in the clinical setting, we conducted drug prediction analyses of noteworthy targets. The primary aim of this study was to identify promising protein biomarkers and genes that could be causally linked and could serve as drug targets for the development of future treatments for DN. This study offers novel insights and methodological approaches that contribute to advances in precision medicine for DN management.

MATERIALS AND METHODS

Study design

As shown in Figure 1, we implemented a seven-step analytical pipeline.

Figure 1 Overview of the study design in our Mendelian randomization and colocalization study. DN: Diabetic neuropathy; UKB-PPP: United Kingdom biobank pharma proteomics project; eQTL: Expression quantitative trait loci; GTEx: Genotype-tissue expression; GWAS: Genome-wide association studies; MR: Mendelian randomization; SMR: Summary-data-based Mendelian randomization; FDR: False discovery rate; GO: Gene Ontology; Cis-pQTL: Cis-protein quantitative trait loci; HEIDI: Heterogeneity dependency tool; PheWAS: Phenome-Wide Association Study.

Proteome-wide target discovery: Candidate protein targets for DN were identified using proteome-wide association studies, based on two-sample MR with protein quantitative trait loci (pQTL) instruments derived from large genome-wide association datasets.

Functional context: Potential biological functions of these proteins were characterized through protein-protein interaction (PPI) network analysis and Gene Ontology enrichment.

Transcript-level validation: Expression of genes encoding the candidate plasma proteins were evaluated using summary data-based MR (SMR), thereby assessing the consistency between protein and transcript associations.

Locus sharing: A co-localization analysis was performed to test whether DN and the protein-encoding genes share the same causal variants.

Evidence grading: MR, SMR, and co-localization findings were integrated to assess evidential strength and prioritize Tier 1 targets.

Therapeutic tractability and safety: For Tier 1 targets, candidate drug prediction and molecular docking were conducted to evaluate therapeutic potential, and phenome-wide association studies were used to explore possible target side effects.

Mechanistic mediation: A two-step mediation MR framework was applied to investigate whether immune cell subsets mediated the effects of Tier 1 targets on DN.

Data source of plasma proteome

The pQTL linked to plasma proteins were derived from two of the most extensive independent GWAS (Figure 1).

deCODE: Utilizing the SomaScan platform, we identified 28191 genetic associations (P < 1.8 × 10^-9) for 4907 aptamers in a cohort of 35559 Icelanders. The data were sourced from two primary projects: The Icelandic Cancer Project (52% of participants) and various genetic initiatives at deCODE Genetics, Reykjavík, Iceland (48% of the participants). Utilizing precalculated summary statistics, recursive conditional analysis was used to identify the most impactful variant within each region (± 1 Mb) as the primary pQTL (n = 18084), while other variants were categorized as secondary pQTLs (n = 10107)[11].

UKB-Pharma Proteomics Project: 23588 primary (sentinel) genetic associations (P < 1.7 × 10^-11, clumping ± 1 Mb, r² < 0.8) were identified for 2923 proteins in 54219 participants from the UKB Pharma Proteomics Project (PPP) utilizing the Olink platform[12]. This GWAS confirmed that 84% of the pQTLs identified in prior antibody-based research replicated 38% of the pQTLs identified in aptamer-based studies[12].

For each plasma protein in the two sets of pQTL data, single nucleotide polymorphisms (SNPs) with minor allele frequencies of at least 1% and genome-wide significance (P < 5.0 × 10^-8) were preserved. Furthermore, SNPs in high linkage disequilibrium (LD) (r² > 0.01 in the 1000 Genomes Project) were considered redundant and omitted from the analysis.

Data sources of transcriptome

Gene expression data were obtained using the GTEx project (version 8) comprising 15201 samples. Plasma- and tissue-specific cis-expression quantitative trait loci (eQTLs) from GTEx (version 8) were used to explore tissue-specific correlations and the potential unintended effects of medications targeting specific genes[20].

Gene expression was validated using the eQTLGen Consortium, which offers a sizable sample size of 31684 individuals to discover SNPs linked to the expression of genes associated with plasma protein targets[19]. This extensive dataset facilitated the validation of the findings from GTEx (version 8).

Data sources of DN

The primary data were sourced from the ninth release of the FinnGen Consortium. The DN dataset comprises 274660 adult Finnish participants, with 2843 cases and 271817 controls[5]. Additionally, GWAS summary data from another DN investigation conducted by Zorina-Lichtenwalter et al[21] were used. This dataset is accessible in the GWAS catalog and encompasses 435971 adult finish participants (772 cases and 435199 controls). GWAS data from both studies were derived from two independent and non-overlapping samples of European descent (Supplementary Table 1).

Data sources of 731 immune cells

This study analyzed 731 immune phenotypes, including 118 absolute cell counts, 389 median fluorescence intensities, 32 morphological parameters, and 192 relative cell counts. Data from 3757 European descendants identified 122 important independent associated signals at 70 loci, 53 new loci, and the regulatory mechanisms of 459 cellular features. Four Illumina arrays were used for genotyping, covering approximately 22 million SNPs, which were substituted based on the Sardinia sequence reference panel. Covariate-adjusted association analysis was performed, considering sex, age, and other relevant factors (Supplementary Table 1)[22].

Proteome wide MR analysis

Initially, the extensive plasma proteome dataset from deCODE was screened using a p-value threshold of < 5^e-8. Before conducting the MR analysis, the deCODE protein and DN data were harmonized to remove SNPs associated with DN. Subsequently, we performed an MR Pleiotropy RESidual Sum and Outlier test, which is recommended when the proportion of horizontal pleiotropy is < 50%[23].

Through MR analysis, potential protein targets of DN were identified in a wide range of plasma proteomic datasets. In the preliminary analysis, only the cis-pQTL was used as the instrumental variable (IVs) for each protein, and these results were based on the DN data from individuals of European descent. For proteins linked to a single cis-QTL, the Wald ratio and delta method were applied to estimate odds ratios (ORs) and corresponding confidence intervals (CIs)[24]. For proteins linked to multiple cis-pQTLs, estimates were derived using the inverse variance-weighted method[25].

This method must follow three key assumptions, as shown in Figure 1: (1) The IV is related to exposure; (2) The IV only affects the results through exposure (represented by the red cross); and (3) The IV is independent of confounding factors (also represented by the red cross). We implemented various measures to test the hypotheses. The MR-Egger intercept test was used to evaluate the level of pleiotropy[26]. To strengthen the control over potential reverse causality and pleiotropy, we used Steiger filters to eliminate pQTL that explained more variability than the corresponding proteins in the DN. We can minimize the genetic confusion caused by horizontal pleiotropy and LD by limiting IVs to cis-pQTL and combining it with colocalization analysis[14,27]. Cochran’s Q test was used to examine the heterogeneity. To address the challenges related to multiple testing, significance was determined based on false discovery rate (FDR) adjusted P values.

Reverse causality detection

Steiger filtering was implemented to verify the direction of the relationship between the proteins and the DN[28].

PPI network construction

Evaluation and analysis of the PPI network were conducted to enhance our understanding of protein interactions. PPI network analysis was performed for proteins showing significant associations with DN, with FDR < 0.05. The PPI network was constructed using a confidence score threshold of 0.4, which was the minimum required interaction score[29]. The results of the PPI network analysis were visualized using Cytoscape version 3.9.1[30].

To deepen our understanding of gene function, we conducted gene function annotation of the identified protein-coding genes to elucidate the biological relevance of the candidate genes. Enrichment was considered statistically significant when the Q value was < 0.05.

Transcriptome wide MR analysis

To validate the identified plasma protein targets, we used the SMR method[31] to examine the relationship between the expression of the respective protein-encoding genes in blood samples from GTEx (version 8) and the risk of DN. The SMR method selects the most significantly correlated eQTL SNP as IVs. In addition, the SMR tool combines Heterogeneity Dependency Tool testing to determine whether the observed association between gene expression and outcomes can be attributed to linkage scenarios. The Heterogeneity Dependency Tool test P value < 0.01 indicates a significant association. The primary outcomes were reported as disease ORs per 1-standard deviation change.

Replication and tissue specific analysis

We conducted gene-level investigations using the GTEx blood samples and a tissue-specific dataset of neural tissues for gene expression analyses. The reliability of our results was confirmed by reproducing the gene analysis using blood data from the eQTLGen consortium.

Colocalization analysis

In certain instances, an SNP may reside in regions that influence multiple genes, potentially affecting disease risk through different genetic pathways. Colocalization analysis is instrumental in confirming whether disease risk and pQTL share underlying causal genetic variations. For significant MR findings, we performed colocalization analysis on SNPs located within ± 100 kb of each gene transcription start site associated with DN risk and pQTL using the R software. The set conditions include P1 = 1 × 10^-4, P2 = 1 × 10^-4, P12 = 1 × 10^-4. P1, P2, and P12 represent the likelihood of SNP association with DN, the likelihood of the SNP being an important pQTL, and the likelihood of an association with DN risk and pQTL, respectively. The COLOC software evaluated five hypotheses and quantified the level of support for each hypothesis using posterior probability (PP): PPH0 (unrelated to any trait), PPH1 (related to gene expression but not to DN risk), PPH2 (related to DN risk but not to gene expression), PPH3 (related to DN risk and gene expression but with different causal variations), and PPH4 (related to DN risk and gene expression but with common causal variations). Given the limited power of colocalization analysis, our assessment was restricted to genes with PPH3 + PPH4 ≥ 0.9[16].

Candidate drug prediction

Assessing protein-drug interactions is essential for assessing the potential of target genes as drug targets. The Drug Signatures Database (DSigDB), a comprehensive database, was used[32]. The identified target genes were input into DSigDB, potential candidate drugs were predicted, and the drug activity of the target genes was evaluated.

Molecular docking

To analyze the binding affinity and interaction modes between candidate drugs and targets, we employed a molecular docking platform hosted at https://www.home-for-researchers.com/. Two-dimensional structures of valproic acid (compound CID: 3121) and cyclosporine A (compound CID: 5284373) were retrieved from PubChem[33]. The three-dimensional coordinates of BTN3A1 (Protein Data Bank ID: 4F80) and MICB (Protein Data Bank ID: 1JE6) were selected from the platform. All docking tasks were performed with the default parameters provided by the web server; the complete docking logs, including receptor-ligand poses, grid maps, and estimated binding energies, are publicly accessible at the above website. Because the calculated binding energies were used solely for rapid triage, only the non-bonded interaction forces (electrostatic, van der Waals, and hydrogen bond terms) were extracted and recorded for further discussion.

Phenome-wide association analysis

To further assess the pleiotropy and potential off-target effects of the candidate drug targets, we investigated the AstraZeneca PheWAS Portal, a compendium of gene-phenotype associations generated from UKB electronic health records, questionnaire data, and exome/genome-wide continuous traits. All analyses were performed using the default portal parameters[34]. Representative PheWAS Manhattan plots and phenotype association results are displayed in the Supplementary material; for the full, interactive versions please consult the AstraZeneca PheWAS Portal gene-view pages at https://www.azphewas.com/geneView/6319c068-fd59-46d8-85ee-82d82482eb14/MICB/glr/binary (MICB) and https://www.azphewas.com/geneView/6319c068-fd59-46d8-85ee-82d82482eb14/BTN3A1/glr/binary (BTN3A1).

Mediation analysis of key genes

We conducted a two-step MR study to investigate the association between proteins and the predicted risk of DN and to explore the potential mediation of this relationship by immune cell characteristics. Initially, we evaluated the causal impact of proteins and immune cell characteristics on DN, focusing on those that displayed significant associations with DN. Subsequently, we examined the causal effects of the selected proteins on the specified immune cell characteristics, encompassing the influence of proteins on immune cells (β1), the impact of proteins on DN (β2), and the overall effect of proteins on DN (β3). Estimate the proportion of total effects mediated by each biomarker using the formula: β1 × β2/β3[35].

RESULTS

The relationship between plasma proteins and DN

To strengthen the validity of the MR analysis, we first applied the MR Pleiotropy RESidual Sum and Outlier test to identify and correct for horizontal pleiotropy, which could otherwise bias the causal estimates by violating the IV assumptions (Supplementary Tables 2 and 3). Based on the filtered IVs, MR analysis identified 127 plasma proteins that were nominally associated with DN at a significance threshold of P < 0.05 (Supplementary Table 4). However, after adjusting for multiple testing using the FDR, only nine proteins remained significantly associated with DN at an FDR < 0.05, indicating a more stringent subset of robust causal candidates (Figure 2A). Among these nine proteins, six exhibited negative associations with DN, and three showed positive associations (Figure 2B). Specifically, HSPA1B, PSMB9, BTN3A1, SCGN, NOTUM, and MICB displayed negative correlations with DN, with ORs and 95%CIs of 0.10 (0.00-0.03), 0.04 (0.01-0.20), 0.06 (0.01-0.23), 0.24 (0.12-0.51), 0.54 (0.40-0.74), and 0.69 (0.59-0.81), respectively. Conversely, WARS, BRD2, and CSNK2B were positively linked to DN, with ORs and 95%CIs of 1.56 (1.26-1.93), 1.86 (1.38-2.34), and 4.35 (3.59-5.10), respectively.

Figure 2 Result summary of Mendelian randomization analysis on the associations between plasma proteins and the risk of diabetic neuropathy. A: Volcano plot of Mendelian randomization analysis; B: Forest plot of Mendelian randomization analysis. FDR: False discovery rate.

A sensitivity analysis was performed to assess the robustness of the causal associations. All eight principal target proteins were examined for heterogeneity and pleiotropy, except for MICB, which demonstrated heterogeneity (Table 1), necessitating additional scrutiny to validate its dependability. Furthermore, Steiger filter analyses were used to investigate the reverse causality, ensuring that the identified causality was not affected by the potential reverse effects of DN on these proteins (Supplementary Table 5).

Table 1 Result summary of Mendelian randomization and sensitivity analysis on the associations between plasma proteins and diabetic neuropathy.

Protein	Method	nSNP	Beta	SE	P value	FDR	OR (95%CI)	MR-Egger	Q
								P value	P value IVW
MICB	Weighted median	25	-0.339	0.036	0.000	0.000	0.71 (0.66-0.76)	0.062	0.000
	MR Egger	25	-0.189	0.121	0.132	0.902	0.83 (0.65-1.05)
	Inverse variance weighted	25	-0.371	0.082	0.000	0.002	0.69 (0.59-0.81)
NOTUM	Weighted median	6	-0.526	0.198	0.008	0.473	0.59 (0.40-0.87)	0.126	0.365
	MR Egger	6	-0.067	0.321	0.845	0.993	0.94 (0.50-1.75)
	Inverse variance weighted	6	-0.611	0.161	0.000	0.024	0.54 (0.40-0.74)
WARS	Inverse variance weighted	7	0.445	0.108	0.000	0.008	1.56 (1.26-1.93)	0.781	0.512
	MR Egger	7	0.493	0.198	0.055	0.766	1.64 (1.11-2.42)
	Weighted median	7	0.457	0.123	0.000	0.031	1.58 (1.24-2.01)
HSPA1B	Wald ratio	1	-4.552	0.501	0.000	0.000	0.01 (0.00-0.03)
BTN3A1	Wald ratio	1	-2.852	0.695	0.000	0.009	0.06 (0.01-0.23)
PSMB9	Wald ratio	1	-3.272	0.854	0.000	0.022	0.04 (0.01-0.20)
SCGN	Inverse variance weighted	2	-1.415	0.375	0.000	0.026	0.24 (0.12-0.51)
CSNK2B	Wald ratio	1	1.470	0.888	0.000	0.000	4.35 (3.59-5.10)
BRD2	Wald ratio	1	0.623	0.564	0.000	0.000	1.86 (1.38-2.34)

OR: Odds ratio; CI: Confidence interval; nSNP: Number of single nucleotide polymorphisms; SNP: Single nucleotide polymorphisms; FDR: False discovery rate; MR: Mendelian randomization; IVW: Inverse variance-weighted.

In the replication cohort, observations confirmed significant associations (P < 0.05) between MICB and BRD2 and DN, exhibiting consistent beta values (Table 2). Although CSNK2B and SCGN also passed the P < 0.05 threshold, the directions of their beta values were inconsistent with those in the discovery cohort. Subsequent analyses utilizing the UKB-PPP plasma protein data and Finngen findings revealed that MICB, WARS, and SCGN plasma proteins were significantly associated (P < 0.05) with consistent beta values (Table 3). When utilizing UKB-PPP as exposure and the GWAS catalog as outcome, SCGN and MICB met P < 0.05, although SCGN’s beta value direction differed from that observed in the discovery cohort (Table 4, Supplementary Tables 6-8).

Table 2 Mendelian randomization analysis of associations of plasma proteins with diabetic neuropathy in replicated stage (decode vs genome-wide association studies catalog).

Outcome	Exposure	Method	nSNP	Beta	SE	P value
GCST90129436	HSPA1B	Wald ratio	1	-0.997	1.117	0.372
GCST90129436	PSMB9	Wald ratio	1	1.346	1.542	0.383
GCST90129436	BTN3A1	Wald ratio	1	1.374	1.412	0.331
GCST90129436	SCGN	Inverse variance weighted	2	1.446	0.628	0.021
GCST90129436	NOTUM	Inverse variance weighted	4	-0.056	0.353	0.874
GCST90129436	MICB	Inverse variance weighted	22	-0.188	0.016	0.038
GCST90129436	WARS	Inverse variance weighted	6	-0.033	0.256	0.897
GCST90129436	BRD2	Wald ratio	1	3.393	1.267	0.007
GCST90129436	CSNK2B	Wald ratio	1	-5.066	1.980	0.010

nSNP: Number of single nucleotide polymorphisms; SNP: Single nucleotide polymorphisms.

Table 3 Mendelian randomization analysis of associations of plasma proteins with diabetic neuropathy in replicated stage (United Kingdom biobank pharma proteomics project vs genome-wide association studies Finngen).

Outcome	Exposure	Method	nSNP	Beta	SE	P value
DN	HSPA1B	Inverse variance weighted	2	3.696	2.072	0.074
DN	SCGN	Wald ratio	1	-0.938	0.346	0.007
DN	MICB	Inverse variance weighted	29	-0.187	0.001	0.008
DN	WARS	Inverse variance weighted	3	0.482	0.168	0.004
DN	BRD2	(-)	-	-	-	-
DN	CSNK2B	(-)	-	-	-	-
DN	PSMB9	(-)	-	-	-	-
DN	BTN3A1	(-)	-	-	-	-
DN	NOTUM	(-)	-	-	-	-

(-) represents that there is no match in the dataset. nSNP: Number of single nucleotide polymorphisms; SNP: Single nucleotide polymorphisms; DN: Diabetic neuropathy.

Table 4 Mendelian randomization analysis of associations of plasma proteins with diabetic neuropathy in replicated stage (United Kingdom biobank pharma proteomics project vs genome-wide association studies catalog).

Outcome	Exposure	Method	nSNP	β	SE	P value
GCST90129436	HSPA1B	Inverse variance weighted	3	-	0.61	0.320
				0.60	0
				6
GCST90129436	SCGN	Wald ratio	1	1.322	0.627	0.035
GCST90129436	MICB	Inverse variance weighted	27	-0.132	0.062	0.035
GCST90129436	WARS	Inverse variance weighted	2	-0.211	0.215	0.326
GCST90129436	BRD2	(-)	-	-	-	-
GCST90129436	CSNK2B	(-)	-	-	-	-
GCST90129436	PSMB9	(-)	-	-	-	-
GCST90129436	BTN3A1	(-)	-	-	-	-
GCST90129436	NOTUM	(-)	-	-	-	-

(-) represents that there is no match in the dataset. nSNP: Number of SNP; SNP: Single nucleotide polymorphisms.

In summary, among the nine primary target proteins, SCGN and MICB consistently met the significance threshold of P < 0.05 in all replication cohorts. However, the SCGN displayed some variability in the direction of the beta values, indicating a potentially less stable generalizability. In contrast, the MICB exhibited robust generalizability, highlighting its potential significance as a target.

PPI network analysis and functional enrichment analysis of key target proteins

The nine identified proteins were analyzed using the STRING database to construct the PPI network, which was visualized using Cytoscape. Figure 3A illustrates a 33-node, 445-edge PPI network showing the interactions of these proteins with other proteins. Biological process enrichment analysis, shown in Figure 3B, emphasized two main areas: Inflammation-related processes (such as immune response-activating signal transduction and cytokine-mediated signaling) and neural function-related processes (including neuronal cellular homeostasis and long-term synaptic potentiation). Supplementary Table 9 contains the detailed Gene Ontology terms for these processes, confirming their relevance in the pathophysiological mechanisms of DN.

Figure 3 Results of protein–protein interaction network and Gene Ontology enrichment. A: Protein-protein interaction network built with STRING; B: For Gene Ontology enrichment pathways.

Gene expression associations of key target proteins

By mapping the nine proteins to eight coding genes, we performed SMR analysis using GTEx (version 8) blood data and neurospecific tissues. BTN3A1 and MICB successfully passed the SMR tests (Figure 4A and B), exhibiting consistent associations with a decreased risk of DN across both blood- and neuro-specific tissues (Supplementary Table 10). Validation using eQTLGen blood samples confirmed the directional consistency of the beta values for BTN3A1 and MICB observed in the primary analysis. The comprehensive SMR analysis results are presented in Supplementary Table 10.

Figure 4 Summary data-based Mendelian randomization plot of the MICB and BTN3A1 locus. The top panel depicts diabetic neuropathy genome-wide association studies P values (grey dots), the middle panel depicts expression quantitative trait loci P values, and the bottom panel depicts gene locations on chromosome 6. The red diamonds in the upper panel indicate the summary data-based Mendelian randomization test P values for positive genes, and the dashed lines indicate the significance thresholds. A: A locus plot of MICB; B: A locus plot of BTN3A1. eQTL: Expression quantitative trait loci; GWAS: Genome-wide association studies; SMR: Summary-data-based Mendelian randomization.

Colocalization of key target proteins

Among these proteins, SCGN, HSPA1B, PSMB9, MICB, and BTN3A1 exhibited colocalization with DN (PPH3 + PPH4 > 0.9), whereas CSNK2B, BRD2, NOTUM, and WARS did not meet the specified criteria (Supplementary Table 11).

Integration and comprehensive analysis of key target proteins

To gain a thorough comprehension of the function of key target proteins and to facilitate the exploration of significant targets, we integrated the results of the proteome-wide association studies and TWAS colocalization analyses. Given that proteins represent the ultimate products of gene expression, determining causal connections at the protein level is essential. Consequently, all three levels of candidate causal genes in our classification should provide evidence of a causal association with DN at the protein level. Adhering to this principle, we classified the causal candidate genes into three tiers according to the following criteria: (1) Tier 1 comprises two crucial targets, BTN3A1 and MICB, identified through MR, SMR, and protein colocalization analyses; (2) Tier 2 included genes that passed the colocalization analysis, specifically SCGN, HSPA1B, and PSMB9; and (3) Tier 3 consists of genes that did not pass the colocalization analysis, namely CSNK2B, BRD2, NOTUM, and WARS (Table 5).

Table 5 Summary of key target protein integration results.

Proteins (genes)	Direction for DN risk of protein	Direction for DN risk of gene	Passed SMR analysis	Passed coloc analysis
Tier 1
BTN3A1	(-)	(-)	Yes	Yes
MICB	(-)	(-)	Yes	Yes
Tier 2
PSMB9	(-)	(-)	No	Yes
HSPA1B	(-)	(-)	No	Yes
SCGN	(-)	X	X	Yes
Tier 3
CSNK2B	(+)	(+)	No	No
BRD2	(+)	(+)	No	No
WARS	(+)	X	X	No
NOTUM	(-)	X	X	No

(+) represents that the beta value is positive; (-) represents that the beta value is negative. X represents the loss in the analysis due to data issues, where the relevant gene was not found. “Yes” indicates that the analysis has passed, while “No” indicates that it has not passed. DN: Diabetic neuropathy; SMR: Summary-data-based MR.

Comprehensive analysis of Tier 1 target proteins

Based on this analysis, BTN3A1 and MICB were identified as the prominent targets for DN. To comprehensively investigate their roles, we conducted candidate drug prediction, molecular docking studies, and PheWAS analyses to explore the potential side effects beyond DN.

Here, we identified potentially effective interventional drugs using the DSigDB database. The drugs associated with BTN3A1 and MICB are shown in Supplementary Table 12. The overlap of these drugs revealed that benzo(a)pyrene, cyclosporin A, and valproic acid were significant targets of both BTN3A1 and MICB (Figure 5A and B). However, benzo(a)pyrene, a highly carcinogenic polycyclic aromatic hydrocarbon, is unsuitable for use as a drug and is consequently excluded.

Figure 5 Drug prediction molecular docking. A: Drug prediction molecular docking BTN3A1 and MICB-related drugs; B: Molecular docking plots of BTN3A1 with cyclosporin A and valproic acid; MICB with cyclosporin A and valproic acid; C-F: Predicted binding conformations of cyclosporin A and valproic acid with BTN3A1 and MICB.

To assess the druggability of BTN3A1 and MICB targets, we performed molecular docking investigations to assess the binding affinities of cyclosporin A and valproic acid for these targets. Employing AutoDock Vina, we determined the binding sites and interactions of cyclosporin A and valproic acid with the BTN3A1 and MICB proteins, generating binding energies for each interaction. The results depicted in Figure 5C-F indicated the successful docking of the proteins with the drugs. Every drug candidate formed observable hydrogen bonds and robust electrostatic interactions with the corresponding protein target. Furthermore, binding pockets were effectively occupied by drug candidates. As indicated in Table 6, cyclosporin A displayed the lowest binding energy for both key targets (MICB-cyclosporin A: -11.64 kcal/mol, BTN3A1-cyclosporin A: -10.39 kcal/mol), implying highly stable binding. Hence, CsA may be a promising drug candidate for DN.

Table 6 Binding energy for Tier 1 target proteins with their drugs.

Target	Drug	Binding energy (kcal/mol)
MICB	Cyclosporin A	-11.64
MICB	Valproic acid	-4.27
BTN3A1	Cyclosporin A	-10.39
BTN3A1	Valproic acid	-4.14

To evaluate the potential adverse effects of the identified Tier 1 target proteins on other traits, a PheWAS was performed at the gene level. PheWAS outcomes, which signified the relationship between genetically determined protein expression and specific diseases or traits (Supplementary Figure 1), indicated that Tier 1 target proteins did not show substantial connections with other traits at the gene level (genomic correlation P < 5^e-8). This suggests that the risk of side effects from drugs targeting these proteins, as well as the possibility of horizontal pleiotropy in these genes, is low. These outcomes additionally enhance the reliability of the conclusions of this study.

Mediation analysis of Tier 1 target proteins

Enrichment analysis of the key target proteins revealed their primary involvement in immunity, particularly their close association with T-cell proliferation. Based on these findings, we hypothesized that BTN3A1 and MICB influence the pathophysiology of DN via immune cells. To investigate this hypothesis, we performed two-step network MR analysis to assess the mediating effects of 731 immune cell traits on the relationship between these proteins and DN. We identified proteins and immune cell traits linked to DN. The analysis revealed that BTN3A1 and MICB met the FDR threshold in the protein-to-DN MR analysis (Figure 6A, Supplementary Table 13). In the MR analysis of immune cells in DN, human leukocyte antigen (HLA)-DR++ monocyte% leukocytes, HLA-DR on CD14+ CD16- monocytes, and HLA-DR on CD14+ monocytes passed the FDR test (Figure 6B, Supplementary Table 14).

Figure 6 731 association of immune cells with diabetic neuropathy and mediation analysis. A: Forest plot of proteins with diabetic neuropathy (DN) Mendelian randomization analysis; B: Forest plot of immune cells vs DN Mendelian randomization analysis; C: Forest plot of identified proteins with Mendelian randomization analysis of positive immune cell; D: The effect of MICB on DN may be partially mediated by the human leukocyte antigen-DR (++) monocyte% leukocyte, with a mediating effect of 52.84%. DN: Diabetic neuropathy; WR: Wald ratio; nsnp: Number of single nucleotide polymorphisms; SNP: Single nucleotide polymorphisms; FDR: False discovery rate; HLA-DR: Human leukocyte antigen-DR; IVW: Inverse variance-weighted; OR: Odds ratio; CI: Confidence interval; Pval: P value.

In the subsequent step, we assessed the causal effects of the identified proteins on selected specific immune cell traits (Figure 6C). By integrating these analyses, we determined that MICB protein influences DN outcomes through HLA-DR++ monocytes. The calculated proportion of the mediating effect was 52.84% (Figure 6D, Supplementary Table 15).

DISCUSSION

DN can result in decreased QOL among affected individuals[1,36,37]. The primary neural alterations associated with DN include the following: (1) Axonal degeneration: Axonal atrophy, degeneration, and rupture hinder nerve signal transmission, resulting in diminished sensation, numbness, and pain[38,39]; (2) Demyelination: Injury to the myelin sheath diminishes its protective and insulating properties, slowing nerve conduction and disrupting action potential transmission[40,41]; (3) Microvascular changes: Diabetes-induced damage to the microvasculature affects the nerve blood supply, leading to inadequate nourishment and worsening nerve damage[42-45]; and (4) Metabolic dysregulation: Aberrant metabolic processes within nerve cells, such as impaired energy metabolism, affect nerve cell functionality and viability[46-48]. These alterations collectively manifest as various clinical symptoms, including sensory abnormalities, motor dysfunction, and autonomic dysregulation in individuals with DN.

Prior to DN onset, persistent hyperglycemia in individuals with diabetes has been suggested to trigger immune system dysregulation and dysfunction. Immune dysfunction precedes DN and plays a key role in its progression. The key immune processes implicated in DN include the following: (1) Activation of inflammatory cells: Elevated blood sugar levels induce the activation of inflammatory cells, leading to the release of various inflammatory mediators. These mediators exert direct cytotoxic effects on nerve cells, resulting in nerve damage[49,50]; (2) Autoimmune reactions: Diabetes can initiate autoimmune responses targeting neural components, generating autoantibodies that specifically attack nerve gangliosides and other neural elements, resulting in immune-mediated assaults and subsequent nerve damage[51-53]; (3) Immune cell infiltration: The infiltration of immune cells, including T cells, into nerve tissues affected by DN results in the release of inflammatory mediators, intensifying neural inflammation and contributing to nerve damage[54]; and (4) Oxidative stress and immune interactions: The oxidative stress triggered by hyperglycemia activates the immune system, whereas inflammatory factors stemming from immune responses enhance oxidative stress, establishing a detrimental cycle that exacerbates nerve damage[55-57]. Mitigating oxidative stress, enhancing mitochondrial function, and optimizing responses to nutritional factors can enhance neural function and facilitate DN regeneration[58,59]. Immune responses in DN are intricate and involve interactions among diverse immune cells and molecules, collectively contributing to nerve damage and dysfunction.

Additionally, recent studies have emphasized the pivotal roles of macrophage polarization and inflammasome activation in DN progression. Peripheral nerve-infiltrating macrophages can adopt either a pro-inflammatory M1 or an anti-inflammatory M2 phenotype. Hyperglycemia and advanced glycation end products favor M1 polarization, which increases the release of cytokines such as tumor necrosis factor-α and interleukin-1β, aggravating nerve injury[60-62]. Notably, lncRNA HCG18 promotes M1 polarization through the miR-146a/TRAF6 axis, thereby accelerating diabetic peripheral neuropathy, whereas downregulation of TRAF6 reverses this effect[60]. In contrast, M2 macrophages exert protective effects against DN pain. M2 macrophages dampen inflammation and facilitate regeneration through the secretion of anti-inflammatory mediators and upregulation of nerve repair-related molecules such as p75NTR, ultimately alleviating neuropathic pain[60]. Inflammasome activation is closely associated with the pathogenesis[62]. Salidroside activates adenosine monophosphate-activated protein kinase signaling and suppresses NOD-Like receptor family pyrin domain-containing 3 inflammasome activation, which improves hyperglycemia and insulin resistance, reduces neuroinflammation, and mitigates neuropathic pain in diabetic rats, suggesting that the adenosine monophosphate-activated protein kinase-NOD-Like receptor family pyrin domain-containing 3 axis is a promising therapeutic target[63]. Altogether, strategies that modulate macrophage polarization and inhibit inflammasome activation may limit immune-mediated nerve injury and promote neural repair, thus offering a tractable avenue for DN therapy.

Our study presents these discoveries and affirms the implications of key target proteins in the fundamental biological processes of diabetes, thereby offering a pathological foundation for DN treatment. Using MR analysis, we identified nine significant proteins and ranked their importance, highlighting BTN3A1 and MICB as Tier 1 genes that are likely to be implicated in DN and potential therapeutic targets. BTN3A1 plays a pivotal role in immune system modulation, particularly in T cell functionality[64]. In the context of aneurysm treatment, BTN3A1 has been identified as a promising drug target, with the inhibition of PLAU and PSMA4 potentially reducing aneurysm susceptibility[65]. BTN3A1 is associated with autoimmune conditions, and its gene polymorphisms contribute to disease susceptibility[66,67]. Furthermore, BTN3A1 is a therapeutic target for type 1 diabetes[68]. These observations are in accordance with the outcomes of our study, underscoring BTN3A1’s regulatory role in immune responses and its intimate connection with diabetes, thereby positioning it as a potential target for addressing DN.

MICB, situated in the MHC region, where genes are implicated in various autoimmune disorders, serves a key function in immune surveillance by facilitating the recognition of abnormal cells through its interaction with the NKG2D receptor. Although common variants of MICB have not been directly linked to type 1 diabetes, their genomic location and immune-related activities suggest their potential involvement in other immune-mediated conditions[69]. MICB is speculated to contribute to the onset and progression of DN, underscoring its probable significance as a therapeutic target in light of immune dysregulation observed in DN[14,70].

Although direct evidence linking BTN3A1 and MICB to DN remains limited, our preliminary findings, together with previous literature, suggest that these proteins participate in biological processes that are highly relevant to DN. In other disease settings, BTN3A1 regulates mitochondrial function and modulates the production of reactive oxygen species[71,72]. Because oxidative stress and nerve injury are central to DN pathogenesis[73-75], these functions implicate BTN3A1 as a potential contributor. MICB has been reported to sustain cellular homeostasis by modulating systemic inflammatory responses, stabilizing mitochondrial function, and reducing intracellular oxidative stress[76-79]. Given that mitochondrial dysfunction, oxidative stress, and inflammation are major drivers of DN progression[80-82], BTN3A1 and MICB may emerge as key targets that influence disease development. Building on this rationale, targeted protein degradation (TPD) provides a practical path from genetic prioritization to therapeutic modulation. Proteolysis-targeting chimera (PROTAC) probes, which induce event-driven degradation at catalytic doses, can accelerate target deconvolution for proteins with low abundance or limited ligandability and can validate MR-nominated targets in relevant immune cell subsets. Applying such probes to BTN3A1 and MICB, either directly or by degrading downstream signaling mediators, may clarify the causal mechanisms and improve druggability in DN[83]. Simultaneously, advances that include nano PROTAC formulations, receptor mediated delivery using folate, antibody, or aptamer conjugates, and oligonucleotide guided transcription factor and RNA PROTAC strategies expand the reach of degradation technologies. Lysosome-directed modalities such as lysosome-targeting chimaeras enable the depletion of cell surface proteins that may be pertinent to MICB biology, whereas proteasome-based PROTACs can target intracellular effectors downstream of BTN3A1 and MICB. Editorial perspectives further indicate that TPD applications extend beyond oncology to viral, neurodegenerative, inflammatory, and metabolic conditions, positioning DN as a plausible next area for translation[83,84]. Collectively, these observations support BTN3A1 and MICB as promising therapeutic candidates and justify focused experimental studies to define their roles and mechanisms of action in DN.

Our study delves deeper into the regulatory functions of MICB and BTN3A1 in DN involving the immune cells. Our analysis suggests that MICB confers protective effects against DN through monocytes, offering mechanistic insights into its role in DN. Although our preliminary analyses suggest that MICB and BTN3A1 may play important roles in DN, their cell type-specific expression, subcellular localization, and detailed mechanisms of action remain unclear. In subsequent studies, we will first apply single-cell sequencing to delineate the cell type-specific expression profiles of BTN3A1 and MICB in DN-related tissues. Building on these data, we will manipulate gene expression in selected cell types using CRISPR Cas9 or Cre LoxP systems, as well as targeted overexpression, and then assess the resulting phenotypic changes in disease-relevant models. To define downstream programs, we will conduct high-throughput RNA sequencing and related analyses to identify the regulated target genes and signaling pathways. These studies will provide a more comprehensive understanding of how MICB and BTN3A1 function in DN, including their cell type-specific roles and impact on key signaling cascades, establishing a mechanistic basis for future therapeutic development.

Protein research plays a crucial role in medicine; however, it is not sufficient for understanding its pathological implications. Translating these findings into effective treatments is crucial. The identification of BTN3A1 and MICB as drug targets for DN offers a clear pathway for the development of novel therapies. Using drug prediction and molecular docking, we identified cyclosporin A as a significant target for DN. This immunosuppressant is commonly used to treat inflammatory conditions such as type 1 diabetes and diabetic retinopathy[85-88]. Cyclosporin A has demonstrated efficacy in the early treatment of type 1 diabetes in older children and adults[87,88]. The side effect profile, however, requires further optimization.

Our study provides extensive evidence from various perspectives. This study is the largest and most comprehensive proteomic and transcriptomic MR analysis of drug targets associated with the DN phenotype. Given the rapid progress in TPD delivery methods, expansion of the target landscape, and steps toward clinical translation, integrating PROTAC-based validation and ultimately PROTAC- or lysosome-directed modalities such as lysosome-targeting chimaeras-guided interventions into DN pipelines appears to be feasible and timely for BTN3A1 and MICB.

This study had several limitations. First, the study population was primarily of European descent, which may restrict generalizability. Therefore, validation in other ethnic groups with DN is required. Secondly, most eQTL instruments are derived from whole blood, whereas DN predominantly involves neural tissues, Schwann cells, and dorsal root ganglia. Future studies should integrate tissue-specific eQTL data from human neural tissues with peripheral blood to improve their biological relevance. Third, the limited availability of high-quality pQTLs and eQTLs led to the exclusion of many proteins and coding genes, which constrained both the discovery of additional candidates and the validation of identified targets; alleviating this constraint will require larger sample sizes and broader gene coverage, as well as improved tissue matching in future eQTL resources. Fourth, although appropriate MR methods and sensitivity analyses were employed, residual horizontal pleiotropy could not be ruled out, and MR-based causal signals still required orthogonal experimental confirmation. Fifth, because there is no widely accepted method for comparing the effectiveness of different analytical pipelines, our evaluation of the strength of the evidence relied on consistency across analyses and datasets. Finally, although our results were robust across multiple protein data sources and tissue-specific associations, their translation to clinical practice requires replication in diverse populations and validation through rigorously designed clinical trials.

CONCLUSION

We identified nine DN-associated proteins and prioritized them for DN treatment using MR, SMR, and colocalization analyses. These proteins and their corresponding genes were enriched in pathways related to immunity and nerve tissue, indicating their potential as therapeutic targets for DN. Further experimental and clinical investigations are required to validate these findings and devise effective treatments centered on these targets.