Data collection
The gene expression matrix and the corresponding clinical data for the GSE142025 and GSE96804 datasets, which are related to DN, were retrieved from the Gene Expression Omnibus (GEO) database (https://www.ncbi.nlm.nih.gov/geo/)15. The training set (GSE142025) included kidney tissue samples from 27 patients with DN and nine controls, which were sequenced using the GPL20301 platform. The validation set (GSE96804) comprised sequencing data from 41 DN patients and 20 controls, processed using the GPL17586 platform. Both datasets focus on kidney tissue; the GSE96804 dataset specifically examines the glomerulus, the kidney's primary filtration unit (Figure 1). The GSE142025 dataset (training set) comprises whole kidney tissue samples and provides an extensive overview of the transcriptomic landscape of DN. The GSE96804 dataset (validation set), on the other hand, focuses specifically on glomerular tissue, which is the primary site of glomerular filtration injury. As these two datasets were not directly merged due to differences in platform and tissue, batch effect correction was not applied. Instead, cross-dataset validation was performed independently. Genes with consistent directional changes and an area under the curve (AUC) greater than 0.7 in both datasets were selected as robust candidates, supporting generalizability across kidney compartments.
A total of 1,136 mitochondria-related genes (MRGs) were extracted from the MitoCarta3.0 database (https://www.broadinstitute.org/mitocarta). The active ingredients of LLF were predicted using the Traditional Chinese Medicine Systems Pharmacology (TCMSP) database (http://sm.nwsuaf.edu.cn/lsp/tcmsp.php), based on an oral bioavailability (OB) threshold of ≥30% and a drug-likeness (DL) threshold of ≥0.18. Subsequently, potential target genes for the active constituents were predicted using the Swiss Target Prediction database (http://www.swisstargetprediction.ch/).
Differential expression analysis
Differential expression analysis of GSE142025 (DN vs. control) was performed using the limma package (v3.54.1), with significance criteria set at P.adj < 0.05 and |log2FoldChange| > 0.516 . Volcano plots and heat maps were visualized using the ggplot2 (v 3.3.6) and ComplexHeatmap (v 2.14.0) packages17,18, respectively. DEGs, MRGs, and potential target genes of the active ingredients were intersected, and the overlapping genes were defined as candidate genes. The network linking the active ingredients and candidate genes was constructed using Cytoscape software (v 3.9.0)19.
Functional enrichment analysis and protein-protein interaction (PPI) network construction
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses for candidate genes were performed using the clusterProfiler package (v 4.6.2) to explore their biological functions and associated signaling pathways (P.adjust < 0.05). The candidate genes were then input into the STRING database (https://cn.string-db.org/) to retrieve PPI relationships (confidence level ≥ 0.4), followed by the construction of a PPI network using Cytoscape (v 3.9.0)20.
Machine learning
Four machine-learning algorithms, including random forest (RF), k-nearest neighbor (KNN), partial least squares (PLS), and support vector machine with a radial basis kernel (SVM), were implemented using the caret package (v6.0-93) based on the GSE142025 dataset. The candidate genes identified in the preceding analysis were used as predictor variables, and disease status (DN or control) was used as the outcome. For the KNN model, 10-fold cross-validation was implemented using the trainControl function, with tuneLength = 10. The RF model was fitted with 20 trees (ntree = 20); the PLS and SVM models were fitted using the settings implemented in the original caret workflow.
Model residual distributions were evaluated using the DALEX package (v2.4.3). Permutation-based variable assessment was subsequently performed using the DALEX variable_importance function, with root mean square error (RMSE) specified as the loss function. The resulting dropout_loss represents the RMSE obtained after permutation of the corresponding variable21,22. In the original analytical workflow, variables with a dropout loss < 0.281 across all models were retained, and genes common to all four models were defined as consensus feature genes for subsequent validation.
The four machine-learning models were used primarily for feature prioritization rather than for constructing a final clinical classifier. Accordingly, diagnostic discrimination was subsequently evaluated at the individual-gene level using ROC analysis in both the discovery and validation datasets.
Identification of candidate biomarkers
In GSE142025 and GSE96804, expression differences in feature genes between DN and control samples were evaluated using the Wilcoxon test. Only genes with significant differential expression (P < 0.05) and concordant expression trends across the two datasets were selected for receiver operating characteristic (ROC) curve analysis. The pROC package (v1.18.0) was used to generate ROC curves and calculate the area under the curve (AUC), with genes showing AUC > 0.7 in both datasets classified as candidate biomarkers23.
Gene set enrichment analysis (GSEA)
The biological functions and signaling pathways associated with the candidate biomarkers were further explored using GSEA on the GSE142025 dataset. First, Spearman correlation analysis of candidate biomarkers with all other genes was performed using the psych package (v2.2.9)24. Correlation coefficients were calculated and ranked (from high to low). The reference gene set was c2.cp.kegg.v2023.1.Hs.symbols.gmt from the Molecular Signatures Database (MSigDB, https://www.gsea-msigdb.org/gsea/msigdb/). GSEA was then conducted to evaluate the enrichment of ranked genes in the background gene set using the clusterProfiler package (v4.6.2). Multiple testing correction was applied via the FDR method, and adjusted P-values (denoted as P.adjust) were considered significant if < 0.05.
M6A modification analysis
To investigate RNA methylation modifications of candidate biomarkers, the SRAMP database (http://www.cuilab.cn/sramp/) was utilized to predict m6A modification sites on candidate biomarkers, focusing on high-confidence positions within their secondary structures. The ENCORI database (https://starbase.sysu.edu.cn/) was then used to identify m6A-modified proteins that interact with candidate biomarkers, using the |HepG2 (shRNA)| > 1 parameter to screen key proteins. The RPISeq database (http://pridb.gdcb.iastate.edu/RPISeq/) was subsequently used to predict the likelihood of interactions between key proteins and candidate biomarkers. RNA sequences for both were uploaded in plain text format to generate RF and SVM classifier prediction scores. An interaction was considered significant when the score exceeded 0.525. SRAMP analysis was performed with the ‘High’ threshold for m6A site prediction, using ‘Transcript’ mode and default parameters. ENCORI analysis used the ‘miRNA-mRNA’ interaction function with the ‘HepG2 (shRNA)’ parameter > 1. RPISeq analysis used the RF classifier with default parameters; scores > 0.5 indicated positive interaction. These are computational predictions, not experimental evidence of m6A modification or RNA-protein interactions in kidney tissue. The HepG2 shRNA criterion was derived from ENCORI’s precomputed datasets and may not reflect renal-specific regulation.
Immune infiltration analysis
The CIBERSORT algorithm was applied to estimate the proportions of 22 immune cell types in both control and DN samples from GSE142025, with visualization via a heatmap generated using the ggplot2 package (v3.3.6)26. CIBERSORT was run using the LM22 signature matrix, with 1,000 permutations and quantile normalization disabled (as recommended for microarray data). Samples with CIBERSORT p < 0.05 were retained for further analysis. CIBERSORT estimates immune cell fractions from bulk kidney tissue expression, which cannot resolve compartment-specific infiltration (e.g., glomerular vs. tubulointerstitial) nor distinguish infiltrating leukocytes from resident immune cells. Thus, the reported correlations are at the tissue level and should be validated by spatial methods. Spearman's correlation analysis between differentially abundant immune cells and candidate biomarkers was performed using the psych package.
Construction of networks and molecular docking
MicroRNAs (miRNAs) interacting with candidate biomarkers were predicted using the miRNet database (https://www.mirnet.ca). Subsequently, long noncoding RNAs (lncRNAs) targeting the identified miRNAs were predicted through the TarBase (http://www.diana.pcbi.upenn.edu/tarbase) and starbase (http://starbase.sysu.edu.cn/) databases. lncRNAs common to both databases were selected for network construction. A lncRNA-miRNA-mRNA regulatory network was then built using Cytoscape software. Potential active ingredients targeting candidate biomarkers were chosen to construct an active ingredient-biomarker network. Additionally, active ingredients, candidate biomarkers, and pathways identified in GSEA were incorporated into Cytoscape to create an active ingredient-biomarker-pathway network.
Molecular docking analysis was performed to assess the binding affinity between potential active ingredients and candidate biomarkers. The 3D structures of biomarker proteins were obtained from the Research Collaboratory for Structural Bioinformatics Protein Data Bank (RCSB PDB, https://www.rscb.org/pdb) in PDB file format. The 2D structures of potential active ingredients were retrieved in SDF format from the PubChem database (http://pubchem.ncbi.nlm.nih.gov). Molecular docking was conducted using the CB-Dock platform (http://clab.labshare.cn/cb-dock/php/blinddock.php). A binding energy of less than -5 kcal/mol indicated a strong binding affinity27.
Preparation and authentication of Ligustri Lucidi Fructus
Here, Ligustri Lucidi Fructus (LLF) refers to the dried ripe fruit of Ligustrum lucidum W. T. Aiton (Oleaceae). The botanical material was authenticated according to the Chinese Pharmacopoeia, and a voucher specimen No. 20240506,20240911,20241103 was deposited at Shanxi University of Traditional Chinese Medicine.
For decoction preparation, 200 g of qualified LLF slices were soaked in 1,000 mL of distilled water for 30 min at room temperature. The mixture was boiled vigorously and then simmered gently for 60 min. The filtrate was collected, and the residual herb materials were decocted repeatedly with another 1,000 mL of distilled water for 60 min. The two filtrates were combined, filtered, centrifuged, and concentrated under reduced pressure to obtain a final stock concentration of 1 g of crude drug/mL (total volume of 100 mL). The prepared decoction was stored at 4 °C for short‑term use or −20 °C for long‑term preservation. The quality of LLF and its decoction was strictly identified and verified according to the standards of the Chinese Pharmacopoeia to ensure experimental reliability and reproducibility.
For qualitative identification, thin-layer chromatography was performed. Briefly, an appropriate volume of prepared decoction was centrifuged, and the supernatant was extracted with methanol. After filtration, the sample solution and specnuezhenide reference standard solution were spotted on the same silica gel G plate. After developing, drying, and inspecting under ultraviolet light, the spot of the sample solution showed consistent color and position with the reference compound, confirming the existence of the characteristic bioactive component of LLF.
For quantitative quality control, high‑performance liquid chromatography detection was carried out. The analysis was performed using a C18 column with methanol–water as the mobile phase. The detection wavelength was set at 224 nm. The content of specnuezhenide in the decoction was determined based on the standard curve. The results demonstrated stable and uniform chemical composition of the prepared decoction, ensuring consistent drug quality throughout the animal intervention experiment.
Animal experiments
Twelve SPF-grade male db/db mice (8–9 weeks old) and six age-matched db/m mice were maintained in the SPF animal facility of Shanxi University of Traditional Chinese Medicine. Before the experiments, the animals were acclimated for 7 days under a 12 h light/12 h dark cycle with ad libitum access to food and water. The study was approved by the Ethics Committee of Shanxi University of Traditional Chinese Medicine (Approval No. 2022DW167). Animals exhibiting >20% body weight loss, a moribund condition, or an inability to access food or water were humanely euthanized before the planned study endpoint. At the end of the study, all remaining mice were euthanized by intraperitoneal injection of pentobarbital sodium followed by cervical dislocation in accordance with institutional protocols.
Following the acclimatization period, the establishment of a DN model in db/db mice was confirmed by a tail-vein blood glucose level ≥ 16 mmol/L and microalbuminuria, as indicated by a positive urine microalbumin test strip. Upon successful establishment of the DN model, the db/db mice were randomly divided into two groups (n = 6 per group): the DN model group (DN) and the LLF treatment group (Treatment). Moreover, db/m mice (n = 6) were used as the control group (Control). The dosage was selected based on previous pharmacodynamic studies of LLF in diabetic rats and was converted to the human equivalent dose using body surface area normalization28. The Control and DN groups were given distilled water, whereas the treatment group received 3.5 g/kg LLF for 8 weeks. After 8 weeks of administration, all mice were euthanized to collect serum, urine, and kidney tissue for subsequent testing.
Blood and urine indicators
Serum glucose levels were analyzed using a fully Automatic Blood Biochemistry Analyzer. The urinary microalbumin concentration was measured according to the kit instructions (Supplemental File 1).
Pathological observation of mouse kidney tissues
Kidney tissues were processed for histopathological examination. After fixation in 4% paraformaldehyde, the tissues were washed, dehydrated, embedded in paraffin, and sectioned. Hematoxylin and eosin (HE) staining was then performed, and the stained sections were examined under an optical microscope to evaluate pathological changes.
Reverse transcription quantitative polymerase chain reaction (RT-qPCR)
RT-qPCR was used to determine the expression of the candidate biomarkers in mouse kidney tissue. Total RNA was extracted according to the manufacturer's instructions, after which RNA concentration and quality were assessed (Table 1). cDNA was synthesized from the extracted RNA using the cDNA Synthesis Kit. Amplification was performed with the primer pairs listed in Table 1, using GAPDH as the reference gene. Relative expression levels were calculated with the 2−ΔΔCt method13,26.
Statistical analysis
All statistical analyses were carried out using R software (version 4.2.2) together with the software required for the corresponding experimental procedures. Unless otherwise indicated, all statistical tests were two-sided, and differences were considered statistically significant at P < 0.05. Transcriptomic differential-expression analysis was performed using the limma package. Genes with an adjusted P < 0.05 and an absolute log2fold change greater than 0.5 were defined as differentially expressed.
For comparisons of candidate-gene expression between independent DN and control samples, nonparametric Wilcoxon rank-sum tests were used where indicated in the original analytical workflow. Correlations between candidate biomarkers and immune cell fractions were assessed using Spearman's rank correlation.
Experimental data are presented as mean ± SD. Comparisons among the three independent animal groups were performed using one-way analysis of variance when the assumptions of parametric analysis were met. For post hoc comparisons, the least significant difference test was used when variances were homogeneous, whereas Dunnett's T3 test was used when variances were unequal. The baseline and week-8 pharmacodynamic measurements were analyzed and presented separately; no inference regarding a group-by-time interaction was made. The R scripts and source data used for bioinformatics and machine-learning analyses are provided in Supplemental File 1.