Datasets of gene expression studies focusing on MCD were selected from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). These datasets were required to explicitly include renal tissue samples from MCD patients and healthy controls with complete and downloadable data, and their sample sizes had to meet the requirements for statistical analysis (training set sample size ≥20; validation set sample size ≥15). Ultimately, three datasets were selected: GSE216841, GSE246204, and GSE139061. The training set GSE216841 contained 35 samples. After excluding 13 samples of idiopathic membranous nephropathy, renal tissue samples from 14 patients with MCD and 8 healthy controls were selected for analysis. Kidney sample data of 12 MCD patients were obtained from GSE246204 (sequencing platform: GPL20301), and the data of 9 healthy control kidneys were obtained from GSE139061 (sequencing platform: GPL20301). The GSE139061 and GSE246204 data were merged into the validation dataset. A total of 3,278 mitochondrial dysfunction-related genes (MDRGs) were retrieved from the GeneCards database (https://www.genecards.org/) (relevance score >5) (Niu et al., 2024) (Supplementary Table S1). The structures of the icariin (ICA) compounds were obtained by searching in the PubChem database (https://pubchem.ncbi.nlm.nih.gov/) using the keyword “ICA,” and the chemical structures of the ICA compounds were uploaded to PharmMapper (http://lilab-ecust.cn/pharmmapper/), TargetNet (http://targetnet.scbdd.com/), the Comparative Toxicogenomics Database (CTD) (https://ctdbase.org//), the Encyclopedia of Traditional Chinese Medicine (ECTM) (http://www.tcmip.cn/ETCM2/front/#/), and the Traditional Chinese Medicines Systems Pharmacology Platform (TCMSP) (https://old.tcmsp-e.com/tcmsp.php) databases to predict ICA targets, remove duplicate targets, and convert target proteins to the genes based on UniProt IDs in the UniProtKB database (https://www.UniProt.org/). The toxicological parameters of ICA were determined using the ProTox II website (https://tox-new.charite.de/protox_II/).

Differentially expressed genes (DEGs) between the MCD and control samples were screened from GSE216841 using the DESeq2 package (v 1.42.0) (Love et al., 2014). Gene names were standardized (probe IDs were converted to official gene symbols via the UniProt database, and duplicate or unannotated probes were removed), and low-expression genes were filtered out (only genes with expression levels >0 in ≥75% of samples were retained) to reduce background noise. The thresholds used were p < 0.05 and |log₂-fold change (FC)| > 0.5. The ggplot2 package (v 3.3.2) (Gustavsson et al., 2022) was used to plot the volcano plot of these DEGs and mark the top ten upregulated and downregulated DEGs. The pheatmap package (version 0.7.7) (Gu et al., 2016) was used to plot a heatmap.

The VennDiagram package (version 1.7.3) (Mao et al., 2022) was used to screen the differential MCD-ICA target genes by intersecting the DEGs, MDRGs, and ICA target genes. The clusterProfiler package (version 3.16.0) (Wu et al., 2021) was used for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses of differential MCD-ICA target genes (p < 0.05).

The interplay of the differential MDRGs we determined was explored, and the differential MDRGs were input into the STRING website for analysis. Cytoscape software (version 3.10.2) (Shannon et al., 2003) was used to visualize the PPI network (https://string-db.org/) (confidence = 0.4). The four algorithms (MCC, stress, MNC, and degree) in the cytoHubba plugin in Cytoscape software were used to screen for candidate genes with the top 25 scores for each algorithm. The genes screened using these four algorithms were then intersected using the VennDiagram package.

Random forest (RF), neural network (NNet), extreme gradient boosting (XGBoost), and support vector machine (SVM) models were built using the caret package (version 6.0–91) (Zhang et al., 2019). The residuals of the RF, NNet, XGBoost, and SVM models were analyzed using the DALEX package (version 2.4.3) (Guan et al., 2023). The classification efficacy of the four models was evaluated based on the area under the curve (AUC) scores in the training and validation sets using the pROC package (version 1.18.0) (Robin et al., 2011). The model with the lowest residuals and the highest AUC was selected as the best prediction model, and the top ten genes ranked using this model were selected as the characterized genes (AUC >0.7).

In the training (GSE216841) and the validation (GSE139061 and GSE246204) sets, the genes with consistent expression trends and significant differences (p < 0.05) between the groups were obtained using the Wilcoxon test for subsequent analysis. The ROC curves of each characterized gene in the training and validation sets were then visualized using the pROC package (version 1.18.0) (Robin et al., 2011), and the characterized genes with AUC values greater than 0.7 were defined as key genes. Correlation analysis of the key genes was performed using the rcorr function in the Hmisc package (version 5.1–3) (http://biostat.mc.vanderbilt.edu/s/Hmisc) (|correlation (cor)| > 0.3, p < 0.05).

A nomogram of the key genes was constructed using the rms package (version 5.1–4) (Xu et al., 2023), and decision curves were visualized using the rmda package (version 1.6) (https://github.com/mdbrown/rmda). Calibration curves were visualized using the regplot package (version 1.1) (Zhang et al., 2018). The objective was to assess the fidelity of the nomogram.

To better understand the biological functions and pathways of the key genes involved in the process of MCD development, “c2.cp.kegg.v2023.1.Hs.symbols.gmt” was obtained from the Molecular Signatures Database (MSigDB) to serve as a background gene set. In GSE216841, the correlation coefficients (p < 0.05) between the key genes and other genes were calculated and ranked using the corrplot package (v 0.92) (Zhang et al., 2023). GSEA was subsequently conducted using the clusterProfiler package (version 3.16.0) (Wang et al., 2022) (adj.p < 0.05, |normalized enrichment score (NES)| > 1). Using the top five KEGG pathways of the key genes, an ICA–key gene–pathway network diagram was constructed using Cytoscape software (version 3.10.2) (Shannon et al., 2003).

The proportions of 22 infiltrating immune cells between the MCD and control groups were analyzed using the CIBERSORT algorithm. The Wilcoxon signed-rank test was used to compare the differences (p < 0.05) in immune-infiltrating cells between the MCD and control groups. Differential immune cells and correlations between key genes and immune cells were determined based on the results of Spearman analysis using the psych package (version 2.2.5) (Robles-Jimenez et al., 2021).

Mitochondrial dynamics genes from the MitoMiner, MitoCarta, and NCBI GEO databases were collected, and the correlations between the 23 mitochondrial kinetic genes and the differential immune cells were determined based on Spearman correlation analysis results conducted using the psych package (version 2.2.5) (Robles-Jimenez et al., 2021).

The Human Protein Atlas (HPA) database (https://www.proteinatlas.org/) was used to analyze the expression of key genes in nine kinds of kidney tissue cells (podocytes, proximal tubular cells, ascending loops of Henle cells, intercalated cells, endothelial cells, fibroblasts, macrophages, T cells, and plasma cells). These were analyzed using the psych package (version 2.2.5) (Robles-Jimenez et al., 2021).

ChIP-X Enrichment Analysis 3 (ChEA3) (https://maayanlab.cloud/chea3/) was used to predict the TFs that targeted the determined key genes. The network depicting the TF–mRNA interactions was established using Cytoscape software (version 3.10.2) (Shannon et al., 2003).

Molecular docking of drugs to the key genes was performed using CB-Dock, and visualization was performed using PyMOL software. The PDB database (https://www.rcsb.org/search/advanced/structure) was searched, and the protein structures of the key genes were downloaded. The molecular structure of ICA was obtained from the PubChem database (https://pubchem.ncbi.nlm.nih.gov/), and the binding energy was <–5 kcal/mol.

Nine frozen tissue samples from rats were collected from the Gansu University of Chinese Medicine, of which three samples were from rats with MCD, three were from control rats, and three were from MCD rats treated with ICA. The ethical approval authority was the Lanzhou Veterinary Research Institute, Chinese Academy of Agricultural Sciences (ethical approval number LVRIAEC-2024–074). Total RNA was extracted from these nine rat tissue samples using TRIzol reagent (Ambion, United States). The primer sequences are detailed in Table 1.

RT-qPCR analysis was conducted on a CFXLFZ006 real-time PCR detection system (Bio-Rad, Shanghai, United States). The collected data were analyzed using the well-established 2^−ΔΔCT method, with GAPDH used as the reference gene for normalization (Livak and Schmittgen, 2001). Finally, GraphPad Prism (version 5.0.0) (Al-Rawi et al., 2023) was used to plot and calculate the p-value.

Kidney tissue samples from SD rats were homogenized, and mitochondria were isolated using differential centrifugation. The homogenate was lysed by adding 100–200 μL of lysis buffer per 20 mg of tissue, followed by conducting thorough homogenization using a glass homogenizer or an equivalent device. Complete tissue lysis was ensured through adequate homogenization. After lysis, the samples were centrifuged at 12,000 × g for 5 min at 4°C, and the supernatant was collected. An ATP assay kit (Beyotime, Cat. #S0026, Shanghai, China) was used according to the manufacturer’s instructions, and the reaction mixture was incubated in the dark at room temperature for 5 min. Absorbance at 560 nm was then measured using a microplate reader. The ATP concentrations in the samples were calculated based on an ATP standard curve.

Kidney tissues were fixed in 4% glutaraldehyde for 4 h, rinsed with phosphate buffer, and post-fixed in 1% osmium tetroxide for 2 h. The samples were then dehydrated through a graded ethanol series, infiltrated with an acetone–epoxy resin mixture, and embedded in pure epoxy resin. Ultrathin sections were obtained from the samples and were then double-stained with uranyl acetate and lead citrate for electron microscopy observation.

Bioinformatics analyses were performed using the R programming language (version 4.2.2). Differences between two groups were determined using the Wilcoxon rank-sum test, and differences between the PCR experimental groups were determined using the Mann–Whitney U test (p < 0.05). The analysis process of this study is detailed in Figure 1.

Data from GSE216841 were analyzed to identify the genes that were differentially expressed between minimal change disease (MCD) and normal samples. A total of 2,297 differentially expressed genes (DEGs) were identified; there were 1,023 upregulated and 1,274 downregulated genes in the MCD samples, with the |log₂FC| values of the top ten upregulated genes (such as NOD2, FOXM1, and CD53) and downregulated DEGs (such as ACTA1, PRTG, and KIT) revealed in the volcano plot and heatmaps (Figures 2A,B).

Next, 797 ICA target genes were generated from the five databases (Supplementary Table S2). The ICA–ICA target gene network map showed the associations between ICA and the 797 predicted target genes. The pharmacological properties, toxicological reports, and 2D and 3D structures of ICA are presented in Table 2. The intersection of 2,297 DEGs, 3,278 MDRGs, and 797 ICA target genes was then used to obtain 45 differential MCD-ICA target genes for subsequent studies (Figure 2C).

In order to understand the function of the 45 differential MCD-ICA target genes, a Gene Ontology (GO) analysis was performed, which revealed that 19 candidate genes were associated with 552 GO signaling pathways, including 487 biological processes, 21 cellular components, and 44 molecular functions, demonstrating significant enrichment of five GO terms, such as the development of striated muscle tissue and the maturation of cardiac muscle tissue and the platelet alpha granule lumen (Figure 2D).

Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed 21 pathways, including the MAPK signaling pathway, the hematopoietic cell lineage, and the arginine biosynthesis (Figure 2E) (p < 0.05).

In order to obtain the candidate genes, 45 differential MCD-ICA target genes were used for protein–protein interaction (PPI) network construction (confidence = 0.4), and the genes with the top 25 scores in each algorithm were screened using four algorithms (MCC, stress, MNC, and degree); the results are shown in Figures 3A–D. The intersection of the genes screened using the four algorithms yielded a total of 19 candidate genes (Figure 3E).

Using the 19 candidate genes, random forest (RF), neural network (NNet), extreme gradient boosting (XGBoost), and support vector machine (SVM) models were built with the training set. A comparison of the cumulative residual distribution plots (Figure 4A) and residual box line plots (Figure 4B) of the four models revealed that the NNet model corresponded to the smallest residual value.

Furthermore, the receiver operating characteristic (ROC) curves of the four models were plotted using the training (Figure 4C) and validation sets (Figure 4D). The area under the curve (AUC) value of the NNet model was greater than 0.7, indicating that the model had good predictive performance. Finally, the top ten genes in the NNet model were selected as feature genes (KIT, XDH, ITGA2B, CSF1R, PLG, ALB, NOS3, IGF1, CTSB, and ANPEP) (Figure 4E).

In order to obtain key genes of diagnostic significance in MCD, three genes (ALB, ANPEP, and XDH) with notable differences in expression and consistent trends between the MCD and control groups were obtained using the Wilcoxon test; these three genes were revealed to be downregulated in MCD (p < 0.05) (Figures 4F,G).

Similarly, reverse transcription quantitative polymerase chain reaction (RT-qPCR) analysis revealed that the expression of ALB, ANPEP, and XDH was obviously lower in the MCD samples than in the control samples and that the expression of ALB and XDH was significantly increased in the MCD samples after ICA intervention (p < 0.05) (Figure 4H).

Next, the ROC curves of the three genes were plotted, and the AUC values of the two genes in the training (GSE216841) and validation (GSE139061 and GSE246204) sets were found to be greater than 0.7. These two genes were defined as the key genes (ANPEP and XDH) and were used for subsequent analysis (Figures 4I,J).

A nomogram allows the visualization of each predictor and its degree of influence on the outcome event. A nomogram model of the determined key genes was, therefore, constructed in this study to predict the probability of MCD (Figure 5A).

The AUC value was 0.991 (Figure 5B), and the DCA indicated that the nomogram model benefited from higher values than individual key genes (Figure 5C); slopes were close to 1 in the calibration curve (Figure 5D), all of which are indicative of the model’s good predictive effect.

In order to clarify the signaling pathways and biological functions associated with the key genes involved in MCD, gene set enrichment analysis (GSEA) revealed the top five pathways with significant enrichment of two key genes. In the single-gene GO enrichment analysis, the first three pathways of the ANPEP and XDH genes were olfactory receptor activity, sensory perception of smell, and sensory perception of chemical stimulus (Figures 6A,B) (p < 0.05).

In the single-gene KEGG enrichment analysis, the ANPEP gene was specifically enriched in the olfactory signaling pathway, olfactory transduction, and sensory perception. The XDH gene was specifically enriched in the olfactory signaling pathway, the olfactory transduction pathway, and the innate immune system pathway (Figures 6C,D) (p < 0.05).

In addition, to investigate the potential mechanisms through which ICA regulated the key genes, the relationships between key genes and KEGG pathways, which are based on the key genes, ICA, and the top five KEGG pathways of the key genes, the ICA–key gene–TOP5 KEGG (ICA–ANPEP–olfactory signaling pathway) pathway network map relationships (Figure 6E) were constructed. Next, the correlation between the key genes was assessed, the results of which revealed a substantial negative correlation between ANPEP and XDH in the training set MCD samples (cor = −0.55, p-value = 0.043) (Figure 6F).

Immune cells are closely associated with the development of MCD (Ghamdi et al., 2020), and how immune cell infiltration occurs differently between the MCD and control groups was investigated. The infiltration scores of 22 immune cells in the samples from GSE216841 were determined using the CIBERSORT algorithm. The top three cells with the highest percentage of immune cells were resting natural killer (NK) cells (Figure 7A). In addition, eight immune cells, such as resting dendritic cells, eosinophils, and macrophages, were notably different (p < 0.05) between the MCD and control groups (Figure 7B).

Spearman analysis revealed positive correlations between regulatory T cells and resting NK cells (cor = 0.85, p-value = 0.001) and between eosinophils and ANPEP (cor = 0.70, p-value = 0.01), suggesting that these cells may be acting synergistically in certain biological processes. Furthermore, evidence confirmed that resting memory CD4 T cells and regulatory T cells (cor = −0.77, p-value = 0.001), along with monocytes and ANPEP, were the most negatively correlated (cor = −0.71, p-value = 0.01) (Figure 7C).

Mitochondrial dynamics are regulated by fusion and fission proteins, both of which are important for organisms. Therefore, to explore the relationship between mitochondrial dynamics and the MCD immune microenvironment, the relationship between mitochondrial dynamics and the MCD immune microenvironment was explored in this study based on 23 genes related to mitochondrial dynamics. It was revealed that eight genes (DNM1L, MIEF2, MUL1, SLC25A46, STX17, MIGA1, MTCH2, and PLD6) were expressed in the training set samples, and DNM1L and differential immune cells (monocytes) were negatively correlated (cor = −0.63, p-value = 0.05), indicating an antagonistic role of DNM1L and monocytes in disease development (Figure 8A). Next, to determine the expression levels of key genes in nine kinds of renal tissue cells (podocytes, proximal tubular cells, etc.), the magnitudes of expression of key genes in renal tissue cells were analyzed. The results of the analysis revealed that ANPEP was expressed only in the proximal renal tubular cells (Figure 8B), and that XDH was expressed in five types of cells, including endothelial cells, fibroblasts, and macrophages (Figure 8C).

Transcription factors specifically recognize the downstream target genes and build transcriptional complexes, thus playing important roles in regulating various biological processes. Therefore, to understand which TFs regulate key genes during disease development, TFs for the key genes were predicted. The results revealed that 98 TFs regulated ANPEP and 39 regulated XDH. Of these, 12 TFs co-regulated the genes XDH and ANPEP. A TF–key gene regulatory network was then constructed (KLF5–XDH and KLF5–ANPEP) (Figure 9A).

In order to further understand the role of ICA in relation to the key genes and determine the binding affinity of the compound ICA for key genes, the protein structures of the key genes ANPEP and XDH, were subjected to molecular docking with the molecular structure of ICA. The docking results revealed that the binding energy of ICA–ANPEP was −9.4 kcal/mol and that of ICA–XDH was −9.0 kcal/mol, both of which were below −5 kcal/mol. This indicates that the selected ICA has high binding affinity for the key genes ANPEP and XDH (Table 3; Figure 9B).

Colorimetric detection analysis revealed that ATP expression was obviously lower in the MCD samples than in the control samples and that ATP expression was significantly elevated in the MCD samples after ICA intervention (p < 0.01) (Figure 10).

Transmission electron microscopy revealed that the ultrastructure of the mitochondria in the renal tissue from the control group was well preserved, with an intact outer membrane and cristae arranged in an orderly and continuous manner. In contrast, the mitochondria in the MCD group exhibited structural disorganization, including a shrunken and fragmented outer membrane and disordered or even absent cristae; some mitochondria also exhibited vacuolar degeneration. Compared with the MCD group, ICA resulted in a more intact outer membrane in the mitochondria (although some membranes remained slightly shrunken), restored crista integrity, and partially recovered the characteristic filamentous mitochondrial morphology (Figure 11).