This study utilized MDD-related datasets obtained from GEO (http://www.ncbi.nlm.nih.gov/geo/). The GSE201332 dataset, serving as the training set, included whole blood samples from 20 patients with MDD and 20 healthy controls ( Supplementary Table 1 ), while GSE52790, used as the validation set, comprised whole blood samples from 10 patients with MDD and 12 healthy controls. MitoCarta3.0 (https://www.broadinstitute.org/mitocarta) and the HAGR database (https://genomics.senescence.info/) were used to extract 1,136 mitochondria-related genes (MRGs) and 866 aging-related genes (ARGs), respectively. A flowchart of the study is provided in Figure 1 .

Differential expression analysis was performed to identify disease-related genes by screening for differentially expressed genes (DEGs) between the MDD and control groups using the Limma package (v 3.44.3) in GSE201332, applying thresholds of an adjusted p-value < 0.05 and |Log₂FC| > 0.5 (10). The intersection of MRGs, ARGs, and DEGs was then analyzed to pinpoint potential biomarkers. GO and KEGG functional enrichment analyses were conducted using the clusterProfiler package (v 4.0.2) (p-value < 0.05, count > 1) (11). The candidate genes were further analyzed using the STRING database to construct a protein-protein interaction (PPI) network with a threshold of 0.2.

For biomarker screening, LASSO logistic regression was performed on the candidate genes using the glmnet package (v 4.0-2) (12), with family = “binomial” to accommodate the binary outcome variable. To ensure reproducibility, set.seed (30) was used to fix random number generation, and a maximum of 5,000 iterations was allowed to ensure algorithm convergence. Ten-fold cross-validation (nfolds = 10) was employed to evaluate model performance and determine the optimal regularization strength (λ). The λ value with the smallest cross-validation error (lambda.min) was chosen for model fitting, and genes with non-zero coefficients were identified as potential biomarkers. The diagnostic performance of these biomarkers was further evaluated through Receiver Operating Characteristic (ROC) curves in both GSE201332 and GSE52790. The ROC curve plotted the false positive rate on the x-axis and the true positive rate on the y-axis, with the Area Under the Curve (AUC) serving as the quantitative measure. An AUC of 0.5 indicated random guessing, an AUC > 0.7 indicated good discriminative ability, and an AUC > 0.9 suggested excellent diagnostic performance. Genes with an AUC > 0.7 were selected for further analysis. Finally, an ANN model based on the biomarkers was developed using the neuralnet package to further assess their diagnostic performance. The ROC curves for the ANN models were also evaluated in both datasets.

To elucidate the regulatory mechanisms and biological functions of biomarkers, gene set enrichment analysis (GSEA) was performed using the clusterProfiler package (v 4.0.2) (11). GO gene sets were sourced from the org.Hs.eg.db database via the gseGO function, and KEGG gene sets were retrieved using the gseKEGG function with the organism parameter set to “hsa”. First, correlation coefficients between each biomarker and the expression levels of all genes were calculated, and genes were ranked based on these coefficients, from high to low. Enrichment scores were determined using the classic permutation method of GSEA (gene set permutation), with significantly enriched gene sets identified using a threshold of p-value < 0.05. In addition, GeneMANIA (https://genemania.org/) was used to predict the genes and functions associated with the biomarkers.

To explore immune cell-related variations, the study investigated the differences in ssGSEA scores for 29 immune cell types between individuals with neurodegenerative diseases (NDD) and control subjects in the GSE201332 dataset (13). The associations between differentially expressed immune cells and the biomarkers were then analyzed. Moreover, to assess the relationship between biomarkers and immunological factors, various immunomodulators and chemokines were retrieved from the ISIDB database (http://cis.hku.hk/TISIDB/), and their correlations with biomarkers were evaluated using Spearman’s method (|cor| > 0.3, p-value < 0.05).

To examine the molecular regulatory mechanisms of the biomarkers, the NetworkAnalyst platform was used to access the “CORE vertebrates” dataset from the JASPAR database (https://www.networkanalyst.ca/). Transcription factors (TFs) with potential binding affinity to the biomarkers were identified using a motif matching score ≥ 800 and a corresponding p-value ≤ 1e-4. miRNAs related to the biomarkers were then predicted using the miRWalk 3.0 (http://mirwalk.umm.uni-heidelberg.de/) and Starbase v3.0 (http://starbase.sysu.edu.cn/) databases, with Starbase requiring a “Pan-Cancer Conservation” score ≥ 3. The miRNAs identified from both databases were intersected to determine the target miRNAs. Next, lncRNAs corresponding to these miRNAs were predicted using the miRNet 2.0 database (https://www.mirnet.ca/), with the screening condition set as CancerNum > 0 in Starbase v3.0. The lncRNAs identified by both databases were intersected to establish the final set of target lncRNAs. Finally, a TF-gene and lncRNA-miRNA-mRNA regulatory network was constructed using Cytoscape software to visually depict these complex molecular interactions.

To assess the effects of chemotherapeutic drugs on biomarkers, a drug-gene network was constructed using the DSigDB database (https://dsigdb.tanlab.org/DSigDBv1.0/). The structural information for the chemotherapeutic drugs in this network was sourced from PubChem (https://pubchem.ncbi.nlm.nih.gov/). Concurrently, protein sequences and functional data for the key genes were retrieved from the Uniprot database (https://www.uniprot.org/), and their three-dimensional structures were obtained from the PDB database (https://www.rcsb.org/). These protein and drug structures were uploaded to the CB-Dock2 platform (https://cadd.labshare.cn/cb-dock2/php/index.php) for molecular docking analysis. The platform standardized protonation states and automatically identified and defined binding sites using its built-in algorithm. Binding affinity between the proteins and ligands was evaluated based on binding energy, with lower values indicating stronger binding. To validate the stability and reliability of the docking results, additional docking experiments were performed using alternative receptor structures of known binders to assess the binding interactions.

To further elucidate the mechanism of action of drugs and evaluate the stability of drug-biomarker complexes as well as the kinetic characteristics of drug binding, molecular dynamics simulations were conducted using GROMACS 2024.4 software. The simulations followed the AMBER99SB-ILDN force field and utilized the TIP3P water model. A cubic system box was set, ensuring the box edges were 1 nm away from the protein edges, and 0.15 mol/L Na⁺/Cl⁻ ions were added to maintain electrical neutrality. Energy minimization was first performed using the steepest descent method. Subsequently, both a heat bath (NVT, with fixed particles, volume, and temperature) and a pressure bath (NPT, with fixed particles, pressure, and temperature) were applied, employing the V-rescale method for temperature coupling. The reference temperature was set to 300 K, with a time step of 2 femtoseconds, and each simulation phase lasted 100 picoseconds. The final molecular dynamics simulation ran for 20 nanoseconds. To quantify binding characteristics, the root-mean-square deviation (RMSD) of backbone atoms in the protein-ligand complex was calculated to assess conformational stability. The root-mean-square fluctuation (RMSF) of protein backbone atoms was analyzed to observe changes in residue flexibility, and fluctuations in total system energy were monitored to evaluate thermodynamic stability. Additionally, the number of hydrogen bonds and their occupancy between the drug and target were counted to quantify the strength of binding interactions. The distance between the small molecule binding site and the amino acid residues of the protein was measured to evaluate binding stability, interaction mechanisms, and conformational changes.

The position of the biomarker on the chromosome was visualized using the RCircos package (v 1.2.2) (14). Gene sequences for the biomarkers were retrieved from the NCBI database, and subcellular localization was assessed using the mRNALocator database. The relationship between the biomarkers and MDD risk was explored using the CTD database.

To further validate biomarker expression in MDD and control groups, expression levels were analyzed in both the training and validation sets, followed by RT-qPCR validation. Five pairs of whole blood samples were obtained from patients with MDD and healthy controls ( Supplementary Table 2 ) at the People’s Hospital of Xinjiang Uygur Autonomous Region for qRT-PCR analysis.

The study cohort consisted of female patients aged 40 to 50 years, diagnosed with depression, who sought care at the Department of Clinical Psychology, People’s Hospital of Xinjiang Uygur Autonomous Region, in June 2024. A control group of healthy females in the same age range was also included. Inclusion criteria were as follows: (1) Diagnosis: Participants met DSM-5 criteria for MDD, confirmed through structured clinical interviews; (2) Symptom severity: A baseline score ≥18 on the 17-item Hamilton Depression Rating Scale (HAMD-17), indicating moderate-to-severe depression; (3) Age: Adults aged 18–65 years; (4) Treatment status: Participants were not receiving any medication or psychological treatment at the time of enrollment. Exclusion criteria were as follows: (1) Comorbid psychiatric disorders: Axis I disorders (e.g., bipolar disorder, psychosis, primary anxiety disorders) or substance use disorders (within 6 months); (2) High suicide risk: Defined by HAMD item 3 score ≥3, recent suicide attempt, or active ideation with intent; (3) Unstable medical conditions (e.g., neurological disorders, uncontrolled diabetes) or medications influencing mood (e.g., corticosteroids). Participants were required to fast overnight for at least 8 hours before blood collection, which was performed between 8:00 AM and 10:00 AM under controlled temperature and lighting conditions. Blood was drawn by a trained phlebotomist using sterile techniques, and aliquots were stored at -20°C for no more than 2 weeks.

This study was approved by the Clinical Research Ethics Committee of the People’s Hospital of Xinjiang Uygur Autonomous Region (KY2024070801), and all patients provided signed informed consent. To validate biomarker expression, total RNA was extracted from the samples using TRIZOL, according to the manufacturer’s instructions. The first strand of complementary DNA (cDNA) was synthesized from 2 μg of total RNA using the SureScript First Strand cDNA Synthesis Kit (Servicebio, Wuhan, China). RT-qPCR was performed with the 2xUniversal Blue SYBR Green qPCR Master Mix (Servicebio, Wuhan, China). The reaction protocol was as follows: 1 minute at 95°C, followed by 40 cycles of 20 seconds at 95°C, 20 seconds at 55°C, and 30 seconds at 72°C. Primer sequences are listed in Table 1 and were validated for specificity using BLAST. GAPDH was used as the internal reference gene. Gene expression levels were calculated using the 2^-△△Ct method (15). Data analysis and visualization were performed using GraphPad Prism 5 (GraphPad Software Inc., USA).

Bioinformatics analysis was performed using R software. Statistical significance was set at p < 0.05. Due to the small sample size, non-parametric tests (Mann–Whitney U test) were used for group comparisons in PCR experiments to ensure the robustness of the results.

A total of 4,041 DEGs were identified in GSE201332, comprising 2,154 upregulated genes and 1,887 downregulated genes ( Figures 2A, B ). Seven candidate genes were further selected by intersecting DEGs with 1,136 MRGs and 866 ARGs ( Figure 2C ). GO enrichment analysis revealed that these candidate genes were associated with 47 functional categories, including mitochondrial outer membrane, organelle outer membrane, and other relevant terms ( Figure 2D ). KEGG pathway analysis highlighted the involvement of these genes in 10 signaling pathways, such as fatty acid degradation and NOD-like receptor signaling ( Figure 2E ). The PPI network analysis identified key PPIs, including MAVS-BCL2 and NBR1-MMT ( Figure 2F ).

In LASSO regression analysis, the lowest error during cross-validation was achieved with a lambda.min of 0.0115, which led to the selection of four biomarkers: SLC25A5, ALDH2, CPT1C, and IMMT ( Figure 3A ). The ROC curve analysis showed AUC values greater than 0.7 in both GSE201332 and GSE52790, indicating that these biomarkers could effectively differentiate between MDD and control samples ( Figures 3B, C ). The ANN model built using these biomarkers demonstrated excellent diagnostic performance in both the training and validation sets, with AUC values of 1 and 0.95, respectively ( Figures 3D, E ).

GO enrichment analysis indicated that SLC25A5, CPT1C, and IMMT were significantly associated with processes such as cellular protein-containing complex assembly, chromatin organization, and chromosome organization ( Figures 4A–C ). Additionally, CPT1C was linked to the detection of chemical stimuli and sensory perception ( Figure 4D ). In KEGG pathway analysis, SLC25A5, CPT1C, and IMMT were implicated in various biological processes, including ATP-dependent chromatin remodeling, neutrophil extracellular trap formation, and protein processing in the endoplasmic reticulum ( Figures 4E–G ). Moreover, CPT1C was involved in pathways such as cortisol synthesis and secretion, and focal adhesion ( Figure 4H ). The gene-gene interaction (GGI) network revealed additional genes related to biomarkers, including ACSS1 and SLC25A6, which are involved in functions such as organelle outer membrane composition, fatty acid transmembrane transport, and other processes ( Figure 4I ). GO and KEGG pathway enrichment analysis results for the four biomarkers are provided in Supplementary Table 3 .

The ssGSEA algorithm revealed significant differences in the scores of 18 immune cells between the MDD and control groups. For example, activated B cells, activated CD8 T cells, and activated dendritic cells (DCs) displayed reduced expression levels in the MDD group ( Figure 5A ). Most of these immune cell types were positively correlated with one another ( Figure 5B ). Additionally, SLC25A5, ALDH2, and IMMT showed positive correlations with several differential immune cells, including activated B cells, activated CD8 T cells, and activated DCs. In contrast, CPT1C exhibited an inverse relationship with most immune cells, except for Immature B cells and Type 1 T helper cells ( Figure 5C ). Correlation analysis with immune factors demonstrated significant associations between the biomarkers and XCL1, CXCL9, CXCL8, CXCL5, CXCL1, and CCL8 ( Figure 5D ).

A total of 27 TFs were predicted in this study. Bioinformatics analysis suggested that STAT1 may target the promoter regions of ALDH2 and SLC25A5, while NKX3–2 may regulate the transcription of IMMT and SLC25A5 ( Supplementary Figure 1A ). Nine target miRNAs were identified by intersecting the predicted miRNAs from the miRWalk and Starbase databases. Based on this, 79 target lncRNAs were predicted, with the following regulatory pairs: HCP5-hsa-miR-27b-3p-SLC25A5, LINC02535-hsa-miR-30b-5p-ALDH2, among others ( Supplementary Figures 1B–D ). Furthermore, drugs corresponding to ALDH2 and SLC25A5 were predicted in the DSigDB database, including four compounds (acetaldehyde, denatured ethanol, nitroglycerin, disulfiram) and two drugs (clodronic acid, butyric acid) ( Figure 6A ). Molecular docking of these drugs with the biomarkers was performed, with ALDH2 (PDB ID: 1nzw) and nitroglycerin showing the most favorable results, exhibiting a docking energy of -6.4 kcal/mol ( Figure 6B ). To validate the docking results, further analysis revealed that the docking energy between ALDH2 and selective serotonin reuptake inhibitors (SSRIs) was 8.2 kcal/mol ( Figure 6C ). When the 3D structure of ALDH2 was replaced with 1CW3, its binding energy with nitroglycerin was 6.3 kcal/mol ( Figure 6D ), which showed minimal change from the original result, further confirming the accuracy of the docking analysis.

This study investigated the conformational changes and energy stability of ALDH2 upon binding to nitroglycerin using 100 ns molecular dynamics simulations. The results showed that the RMSD value of the ALDH2-nitroglycerin system fluctuated between 0.45 and 0.6 nm, indicating that the protein structure reached dynamic equilibrium between 25 and 100 ns and maintained a stable conformation ( Figure 7A ). RMSF analysis revealed that the flexibility of individual residues ranged from 0.05 to 0.4 nm, reflecting local flexibility while ensuring the overall stability of the binding ( Figure 7B ). Energy monitoring demonstrated that the total system energy remained low with minimal fluctuations, and combined with Gibbs free energy landscape analysis, this further confirmed the thermodynamic stability of the complex ( Figures 7C, D ). Hydrogen bond analysis revealed that nitroglycerin formed 1–2 stable hydrogen bonds with the active site of ALDH2, occasionally increasing to 3–4 bonds, highlighting the significance of non-covalent interactions in maintaining binding stability ( Figure 7E ). Additionally, spatial distance monitoring showed that the distances between key binding sites (Residues 150/179) and nitroglycerin stabilized within the ranges of 0.6-0.8 nm and 0.45-0.65 nm, respectively, without a consistent directional change, further confirming the sustained stability of the binding state ( Figure 7F ). In summary, the ALDH2-nitroglycerin complex exhibited stable conformation, favorable thermodynamic properties, and sustained interactions, demonstrating the robustness of their binding.

Chromosomal localization analysis indicated that SLC25A5 is located on chromosome X, ALDH2 on chromosome 12, CPT1C on chromosome 19, and IMMT on chromosome 2 ( Figure 8A ). In subcellular localization analysis, IMMT was localized to the nucleus, while SLC25A5, ALDH2, and CPT1C were localized to the cytoplasm ( Figure 8B ). Furthermore, based on the CTD database, the biomarkers displayed higher scores in depression-related diseases, suggesting that they play a pivotal role in depression pathogenesis ( Figure 8C ).

The expression levels of the biomarkers were further validated. In both the training and validation sets, SLC25A5 and IMMT expression were significantly decreased in the MDD group, while CPT1C showed an opposite expression trend. ALDH2 was downregulated in the MDD group in both datasets, although this change was not statistically significant in the validation set ( Figures 9A, B ). RT-qPCR validation revealed that in patients with MDD, IMMT and SLC25A5 expression were significantly lower, while CPT1C expression was markedly elevated compared to controls (p < 0.05) ( Figure 9C ). These findings were consistent with the bioinformatics analysis. Only ALDH2 expression did not show a significant difference (p > 0.05) ( Figure 9C ).