Methods

Data sources

Gene expression data used in this study were downloaded from the Gene Expression Omnibus (GEO) database. Two datasets were included: the glomerular microarray dataset GSE96804, containing 41 biopsy-confirmed diabetic nephropathy (DN) samples and 20 healthy controls, and the tubular transcriptomic dataset GSE294519, consisting of 23 biopsy-confirmed DN samples and 13 healthy controls. All patients in both datasets were diagnosed by renal biopsy. Differential expression analysis was performed independently for each dataset, after which the DEGs identified from the two renal compartments were integrated to construct a comprehensive DN-related gene set. To construct the FUNDC1 interaction network, key interacting genes were retrieved from BioGRID, STRING, and GeneMANIA using default parameters. Mitochondrial dysfunction–related genes (MDRGs) were obtained from the GeneCards database by querying the keywords “mitochondrial dysfunction” and “mitophagy,” with a relevance score >5 as the selection criterion.

To evaluate the genetic relevance of candidate genes, cis-eQTL data of whole blood were obtained from eQTLGen (n = 31,684, predominantly of European ancestry), and DN GWAS summary statistics were retrieved from the GWAS Catalog (ID: GCST90475670; 25,234 DN cases and 413,872 controls). SMR analyses were performed using these datasets. All summary-level data were derived from European populations, and exposure and outcome datasets originated from independent GWAS consortia to minimize sample overlap and population heterogeneity.

Differential expression analysis

Differentially expressed genes (DEGs) were independently identified from the glomerular dataset GSE96804 and the tubular dataset GSE294519. For GSE96804, when multiple probes corresponded to the same gene symbol, the probe with the highest average expression was retained. Differential expression analysis was performed using empirical Bayes moderation, and P values were adjusted using the Benjamini–Hochberg method. Genes with an adjusted P value (false discovery rate, FDR) < 0.05 and an absolute log2 fold change (|log2FC|) > 0.5 were considered differentially expressed. Volcano plots and heatmaps were generated for visualization.

Identification of differential interaction genes and enrichment analysis

To identify FUNDC1-associated dysregulated modules in DN, DEGs, MDRGs, and FUNDC1-interacting genes were intersected to obtain FUNDC1–DN hub genes. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were subsequently performed.

SMR identification of key genes and functional characterization of candidate genes

Cis-eQTL variants with P < 5 × 10−8 were selected as instrumental variables and harmonized with DN GWAS summary statistics. Linkage disequilibrium pruning was performed using the 1000 Genomes Project reference panel (r2 < 0.001 within a 1000-kb window), and weak instruments (F-statistic < 10) were excluded. SMR analysis was conducted to evaluate the association between genetically predicted gene expression and DN. The HEIDI (heterogeneity in dependent instruments) test was subsequently performed to distinguish pleiotropy from linkage disequilibrium, and genes with SMR P < 0.05 and HEIDI P > 0.05 were considered candidate genes. For significant genes, regional association plots integrating cis-eQTL, GWAS, and SMR results were generated to visualize the local association signals.

For functional exploration, KEGG gene sets (c2.cp.kegg.v2026.1.Hs.symbols.gmt) were downloaded from MSigDB. Spearman correlation coefficients between each key gene and all other genes were calculated in the DN dataset, and genes were ranked according to the correlation coefficients. Gene set enrichment analysis (GSEA) was performed using the ranked gene list, with adjusted P < 0.05 and |NES| > 1 considered statistically significant. Based on the top enriched KEGG pathways, a FUNDC1–hub gene–pathway interaction network was constructed and visualized using Cytoscape. Correlation analysis between key genes was performed using |cor| > 0.3 and P < 0.05 as the significance criteria.

Construction of the transcriptional regulatory network and cross-compartment transcriptomic assessment

Potential upstream transcription factors (TFs) associated with key genes were predicted using ChEA3, which integrates multiple ChIP-seq–based resources. TFs were prioritized based on integrated enrichment scores. Top-ranked TFs were used to construct a TF–mRNA regulatory network and visualized.

To further evaluate the expression patterns of the final candidate genes across different renal compartments, the expression of the identified hub genes was additionally assessed in an independent transcriptomic dataset derived from a complementary renal compartment. The same differential expression analysis pipeline and statistical criteria described above were applied to evaluate the consistency of expression direction and statistical significance across renal tissue compartments.

Statistical analysis and software

All analyses were performed using R (version 4.2.3), SMR (version 1.3.1), and Cytoscape (version 3.10.3). Differential expression analysis of the microarray dataset was performed using the limma package. Data visualization was conducted using the ggplot2 and pheatmap packages. Functional enrichment analyses were performed using the clusterProfiler package, and correlation analyses were conducted using the Hmisc package. Multiple testing correction was performed using the Benjamini–Hochberg method.