Key Aging-Related Genes Linked to Immune Problems in Gum Disease Found Using Gene Analysis and Machine Learning

Gene expression data for the GSE10334 and GSE223924 datasets associated with PD were retrieved from the Gene Expression Omnibus database. Specifically, the GSE10334 dataset (training set) (platform: GPL570) comprised 183 PD patient samples and 64 healthy gingival tissue controls (H. [15]). The GSE223924 dataset (validation set) (platform: GPL24676) included 10 PD patient samples and 10 healthy peripheral gingival tissue controls [16]. Normalization of the GSE10334 dataset was performed by first determining whether the data required a log2 transformation, followed by conventional array normalization at the probe level. After that, the probes were mapped to genes, and duplicates were removed based on average expression, retaining only protein-coding genes. For the GSE223924 dataset, the prenormalized FPKM matrix provided by the GEO database was directly used. Additionally, 81 SASPRGs were obtained from one study [17], and 83 SASPRGs were obtained from another [18]. After data integration and removal of duplicates, 153 SASPRGs were identified in total (Table S1).

To identify differentially expressed genes (DEGs) between PD and control samples in the GSE10334 dataset, differential expression analysis was conducted using the limma package (v 3.58.1) [19] under the thresholds of p adjust < 0.05 and |log2FoldChange (FC)| > 0.5 (D.-B. [20]). Visualization of DEGs was performed with volcano and heat map plots generated by the ggplot2 package (v 3.4.1) [21].

Intersections between DEGs and SASPRGs (candidate genes) were analyzed using the ggvenn package (v 0.1.10) [22]. Gene Ontology (GO) enrichment analysis was employed to functionally annotate genes based on the GO database, encompassing three categories: biological process (BP), cellular component (CC), and molecular function (MF). Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed using the KEGG database to annotate pathways of all identified protein sets or screened differentially expressed proteins, thereby determining the principal metabolic and signal transduction pathways associated with these genes or proteins. GO and KEGG enrichment analyses of candidate genes were conducted with the clusterProfiler package (v 4.8.3) (T. [23]) (p < 0.05), and the results were visualized using ggplot2 (v 3.4.1). Candidate genes were subsequently submitted to the STRING database to construct a PPI network (confidence score ≥ 0.4)[24], and the network was visualized using Cytoscape software (v 3.10.2) (P. [25]). Furthermore, the CytoHubba plugin of Cytoscape software (v 3.10.2) was applied with the top five algorithms, maximal clique centrality (MCC), maximum neighborhood component (MNC), density of maximum neighborhood component (DMNC), degree, and edge percolated component (EPC), to evaluate the significance of each candidate gene. The intersection of the Top 15 proteins derived from each algorithm was designated as the final set of candidate genes for subsequent analysis.

Machine learning, a widely applied analytical approach for gene feature screening, was employed to refine the final candidate genes obtained from the PPI analysis. Three algorithms, least absolute shrinkage and selection operator (LASSO) regression, support vector machine–recursive feature elimination (SVM-RFE), and eXtreme gradient boosting (XGBoost), were applied. LASSO regression (glmnet package, v 4.1-8) (Y. [26]) constructs a penalty function that reduces the coefficients of less informative variables to zero, thereby excluding them from the model. And a 10-fold cross-validation was conducted. Genes exerting greater influence on PD were retained as LASSO genes, and the corresponding regression coefficients were recorded. Based on all samples in the GSE10334 dataset, SVM-RFE (caret package, v 6.0-94) (Z. [27]) was conducted with 10-fold cross-validation to rank the importance of each final candidate gene. Error rates and accuracies of iterative combinations were assessed, and the combination with the highest accuracy was designated as optimal. Genes selected at this stage were defined as SVM-RFE genes. Subsequently, XGBoost (xgboost package, v 2.0.3.1) [28] was applied to the final candidate genes, and genes with high importance scores were identified as XGBoost genes. Finally, the ggvenn package (v 0.1.0) [22] was used to compute the intersection of LASSO, SVM-RFE, and XGBoost genes, with the overlapping genes designated as characteristic genes [29].

To assess the expression of characteristic genes and identify key genes in the GSE10334 and GSE223924 datasets, the Wilcoxon test was performed to evaluate expression differences between PD and control samples (p < 0.05). Results were visualized with boxplots generated using the ggplot2 package (v 3.4.1) (Q. [30]). Genes consistently exhibiting significant differential expression across both datasets were shortlisted as candidate key genes.

To determine the discriminatory capacity of candidate key genes, ROC curve analysis was conducted using the pROC package (v 1.18.5) [31] across all samples in both datasets. The area under the curve (AUC) was applied to measure predictive sensitivity and specificity, with values closer to 1 indicating superior accuracy. Genes with AUC > 0.7 in both datasets were defined as key genes.

An ANN is a machine learning model designed to emulate the learning processes of biological neural networks. To assess the predictive capability of key genes for PD, an ANN model was constructed using the neuralnet package (v 1.44.2) (Y. [32]) with all samples from the GSE10334 dataset. A 10-fold hierarchical cross-validation was implemented when constructing the ANN model. Subsequently, the pROC package (v 1.18.5) [31] was applied to perform ROC curve analysis, and an AUC > 0.7 was considered indicative of favorable predictive performance. In addition, the ANN model was validated in the validation set GSE223924 using the same method.

To elucidate the biological functions and pathways associated with key genes in PD, based on all samples from the GSE10334 dataset, the c2.cp.kegg.v7.4.symbols.gmt (v 7.4) file from MSigDB was used as the reference gene set. The psych package (v 2.4.3) [33] was applied to perform the Spearman correlation analysis between each key gene and all remaining genes, generating corresponding correlation coefficients. These coefficients were subsequently ranked, and the clusterProfiler package (v 4.8.3) was employed to conduct GSEA using the ordered results. Significant enrichment was defined as |NES| > 1 and p < 0.05. Visualization of GSEA results was carried out using the GseaVis package (v 0.1.0) (H. [34]).

To evaluate immune cell infiltration in PD, the CIBERSORT algorithm was applied to estimate the abundances of 22 immune cell types between the PD and control groups in the GSE10334 dataset [35]. The Wilcoxon test was subsequently performed to assess differences in infiltration levels of these 22 immune cell types between PD and control samples (p < 0.05), and violin plots were generated using the ggplot2 package (v 3.4.1) to visualize the results. Additionally, associations between key genes and differentially infiltrated immune cells across all samples in the GSE10334 dataset were analyzed using the Spearman correlation in the psych package (v 2.4.3). A threshold of |cor| > 0.3 and p < 0.05 was adopted, and correlation heat maps were constructed with the ggplot2 package (v 3.4.1).

To elucidate the functions of key genes, they were entered into the GeneMANIA database, and an interaction network with functionally similar genes was constructed to identify and characterize their roles.

To determine upstream regulators of the key genes, the TRRUST database was used to predict transcription factors (TFs) binding to these genes, and Cytoscape software (v 3.10.2) was applied to construct the TF-mRNA regulatory network.

To investigate upstream microRNAs (miRNAs) of the key genes in PD, two databases, TargetScan (http://www.targetscan.org↗) and miRWalk, were independently applied to predict mRNA–miRNA interaction pairs. The intersection of miRNAs predicted by both databases was used to identify key miRNAs. A miRNA–mRNA regulatory network was subsequently constructed with Cytoscape software (v 3.10.2) based on these key miRNAs and their target mRNAs.

To identify potential therapeutic agents for PD, the Comparative Toxicogenomics Database (CTD, https://ctdbase.org/↗) was applied to predict drugs or compounds associated with key genes. The Top 20 candidates were selected to construct a drug–key gene regulatory network using Cytoscape software (v 3.10.2). These selected compounds were defined as key drugs for molecular docking. Three-dimensional molecular structures of the drugs were retrieved from the PubChem database. Subsequently, key genes were queried in the UniProt and PDB databases to obtain the highest resolution receptor structures. Molecular docking simulations were carried out using CB-dock2, and the conformation with the lowest docking score was downloaded. PyMOL (v 2.6) was then employed to visualize docking conformations. A lower docking score (binding energy) indicated stronger docking efficacy, and binding energies below −5 kcal/mol were regarded as evidence of favorable affinity between drugs and key gene targets.

To further validate the results of molecular docking, we conducted MD simulations. Specifically, the GROMACS 2024.4 software [36] was used to perform a 20 ns MD simulation to further verify the rationality and reliability of the molecular docking results. The entire process followed the rules of the AMBER99SB-ILDN force field, and the TIP3P water model was used.

Before the experiment began, two stages of pre-equilibration treatments were carried out. In the first stage, the NVT system was used to simulate the equilibration at 300 K for 100 ps to stabilize the system temperature. In the second stage of equilibration, the NPT system was employed to conduct the simulation at 1 bar for 100 ps to stabilize the system pressure.

In the experiment, the parameters and topology files of the protein and small molecule ligand were generated by the "AMBER14SB" force field and the AMNER gaff force field, respectively. In order to avoid the complex long-range interactions caused by the periodic boundary conditions, the periodic boundary conditions were set and optimized, and the simulation system was filled with water molecules. Meanwhile, to make the simulation system electrically neutral, some solvent water molecules were replaced by Na⁺ and Cl⁻ with a concentration of 0.15 mol/L. In addition, the steepest descent method was used to minimize the energy consumption of the entire system.

During the simulation, we calculated the root mean square deviation (RMSD), root mean square fluctuation (RMSF), and energy to reveal the changes in the position of the conformation during the simulation, the flexibility, and motion intensity of the protein amino acids.

To determine whether the expression of key genes in clinical samples corresponded with bioinformatic predictions, gingival tissue samples were collected from 10 healthy controls and 10 patients with PD recruited from the Second Affiliated Hospital of Guizhou University of Traditional Chinese Medicine. All participants signed and completed a form indicating their informed consent, and the ethical approval agency was the Ethics Committee of the Second Affiliated Hospital of Guizhou University of Traditional Chinese Medicine (Ethics Number: KY20241217). GraphPad Prism 10 was employed to generate plots and calculate p values. Total RNA was extracted from the 20 samples using the TRIzol method (Vazyme Biotech Co., Nanjing, China) and subsequently reverse transcribed into complementary DNA (cDNA) with the Hifair III 1st Strand cDNA Synthesis SuperMix for qPCR kit (Yysen Biotechnology, Shanghai, China). RT-qPCR was then performed using 2× Universal Blue SYBR Green qPCR Master Mix (Servicebio, Wuhan, China), with GAPDH serving as the internal reference for quantifying gene expression. Primer sequences are provided in Table 1. Expression levels of biomarkers between the two groups were calculated using the 2^−ΔΔCt algorithm and compared with the t-test (p < 0.05).

To evaluate protein-level expression of key genes in periodontal tissues, immunohistochemical (IHC) staining was conducted on paraffin-embedded sections from 10 PD patients and 10 healthy controls. Tissue sections were deparaffinized, rehydrated, and subjected to antigen retrieval with citrate buffer. Endogenous peroxidase activity was blocked using 3% hydrogen peroxide. After blocking with 5% bovine serum albumin, the sections were incubated overnight at 4°C with primary antibodies against ICAM1 (Cat# 10831-1-AP), CXCL12 (Cat# 17402-1-AP), and MMP3 (Cat# 17873-1-AP), all obtained from Wuhan Sanying Biotechnology Co., Ltd.

On the following day, sections were incubated with horseradish peroxidase-conjugated secondary antibodies, developed using 3,3⁣′-diaminobenzidine, and counterstained with hematoxylin. Stained slides were visualized under a light microscope, and marker expression was quantified with ImageJ software by calculating the integrated optical density (IOD). Statistical comparisons between groups were performed using the t-test, with p < 0.05 considered statistically significant.

Bioinformatic data were analyzed using R software (v 4.3.1). The Wilcoxon test or t-test was applied to evaluate differences between the two sample groups. A p value < 0.05 was regarded as statistically significant.

Differential expression analysis of the GSE10334 dataset identified 1316 DEGs (p adjust < 0.05 and |log2FC| > 0.5). Among these, 657 genes were upregulated, and 659 genes were downregulated in PD samples compared with controls (Figure 1a,b) (Table S2).

The intersection of 1316 DEGs and 153 SASPRGs yielded 30 differentially expressed SASP-related genes in PD (candidate genes) for subsequent analysis (Figure 2a). GO enrichment analysis revealed that these 30 candidate genes were enriched in 945 GO signaling pathways (p < 0.05), including 903 BP pathways (e.g., leukocyte migration and myeloid leukocyte migration) (Table S3), 5 CC pathways (e.g., collagen-containing extracellular matrix and endoplasmic reticulum lumen) (Table S4), and 37 MF pathways (e.g., cytokine activity and cytokine receptor binding) (Table S5). Enriched biological functions were ranked in descending order by count, and the top five included leukocyte migration, collagen-containing extracellular matrix, and cytokine activity (Figure 2b). KEGG enrichment analysis demonstrated that these genes were associated with 53 pathways (p < 0.05) (Table S6). The Top 10 pathways, ranked by count, included cytokine–cytokine receptor interaction and rheumatoid arthritis (Figure 2c).

A PPI network comprising 30 nodes and 254 interactions was then constructed (confidence score ≥ 0.4), indicating relatively strong connectivity among these proteins. Representative interactions included CXCL1 with CXCL12 and MMP9 with MMP3 (Figure 2d).

The Top 15 proteins identified by MCC, MNC, DMNC, degree, and EPC algorithms were intersected, yielding 8 final candidate genes: MMP1, CCL3, ICAM1, CXCL12, CCL5, CXCL2, IL1A, and MMP3 (Figure 2e).

Using the LASSO algorithm, when log(lambda.min) = −4.3153, five genes were identified: ICAM1, CXCL12, CCL5, IL1A, and MMP3, designated as LASSO genes (Figure 3a,b). SVM-RFE analysis indicated the highest predictive accuracy when the optimal number of variables was eight, yielding ICAM1, CCL5, CXCL12, CCL3, IL1A, CXCL2, MMP3, and MMP1, defined as SVM-RFE genes (Figure 3c). Similarly, XGBoost analysis of all GSE10334 samples identified five high-importance genes, ICAM1, CXCL12, CCL5, MMP3, and CXCL2, referred to as XGBoost genes (Figure 3d).

The intersection of LASSO, SVM-RFE, and XGBoost results produced four overlapping genes: ICAM1, CXCL12, CCL5, and MMP3, denoted as characteristic genes (Figure 3e).

Expression of these characteristic genes was validated using the Wilcoxon test in both GSE10334 and GSE223924 datasets (p < 0.05). In the GSE10334 dataset, ICAM1, CXCL12, CCL5, and MMP3 were significantly upregulated in PD samples (Figure 3f). In the GSE223924 dataset, ICAM1, CXCL12, and MMP3 showed significant expression, with expression trends consistent with those observed in GSE10334 (Figure 3g). Therefore, ICAM1, CXCL12, and MMP3 were considered candidate key genes. ROC curve analysis was subsequently performed across all samples from GSE10334 and GSE223924. In GSE10334, AUC values for ICAM1 and CXCL12 exceeded 0.7, while MMP3 remained below this threshold (Figure 3h). In GSE223924, AUC values for ICAM1, CXCL12, and MMP3 all exceeded 0.7 (Figure 3i). Based on these results, ICAM1 and CXCL12 were identified as key genes for subsequent analyses (AUC > 0.7).

To assess the clinical diagnostic value of key genes, an ANN model incorporating CXCL12 and ICAM1 was constructed using all samples from the GSE10334 dataset to predict PD incidence. The results indicated that CXCL12 and ICAM1 possessed strong diagnostic value for PD (Figure 4a). ROC curve analysis further demonstrated an AUC of 0.883 (95% CI = 0.832–0.943) (Figure 4b), reflecting relatively high predictive accuracy. Comparable results were also observed in the GSE223924 dataset (AUC = 0.990, 95% CI = 0.962–1.000) (Figure 4c,d).

GSEA revealed that CXCL12 was enriched in 75 pathways, including olfactory transduction (Figure 5a) (Table S7), while ICAM1 was enriched in 61 pathways, such as ribosome (|NES| > 1, p < 0.05) (Figure 5b) (Table S8). Moreover, both genes were jointly involved in Leishmania infection, cell adhesion molecules (CAMs), hematopoietic cell lineage, and Parkinson's disease pathways. These results suggest that these pathways may contribute to the pathogenesis of PD.

Immune infiltration analysis of 22 immune cell types in the GSE10334 dataset was conducted (Figure 5c). The Wilcoxon test revealed significant differences in 11 cell types, with activated mast cells, neutrophils, and plasma cells displaying higher infiltration levels in the PD group (p < 0.05) (Figure 5d). Additionally, ICAM1 expression showed a significant negative correlation with resting mast cells (cor = −0.44, p = 8.95e − 12) (Figure 5e) (Table S9).

GeneMANIA analysis identified 20 genes functionally associated with CXCL12 and ICAM1, including ITGAL and FGG. These genes interacted through seven distinct modes of action, such as physical interactions and coexpression, and were closely related to seven biological functions, including cytokine activity and chemokine receptor binding. ICAM1 was primarily involved in leukocyte migration (Figure 6a).

A TF-mRNA regulatory network was constructed using the TF prediction database. This network revealed six TFs targeting CXCL12 and 21 targeting ICAM1 (Table S10). Among these, RELA and NFKB1 simultaneously targeted both genes (Figure 6b). Additionally, a miRNA–mRNA regulatory network was generated using Cytoscape software. CXCL12 was predicted to be targeted by 509 miRNAs in the TargetScan database and 984 miRNAs in the miRWalk database, with 190 overlapping miRNAs identified (Table S11). ICAM1 was predicted to be targeted by 628 miRNAs in TargetScan and 183 in miRWalk, with 59 overlapping miRNAs identified (Table S12). The miRNA–mRNA network further revealed that seven miRNAs, including hsa-miR-3972 and hsa-miR-301b-5p, simultaneously targeted both CXCL12 and ICAM1 (Figure 6c).

CTD analysis predicted 342 compounds or drugs targeting CXCL12 (Table S13), including tamoxifen and estradiol. A total of 774 compounds were predicted to target ICAM1, such as glucosamine and acetaminophen (Table S14). Among these, 174 compounds were predicted to simultaneously target both genes. Construction of a drug–key gene regulatory network using the Top 20 compounds revealed that streptozocin, dietary fats, dexamethasone, and tetrachlorodibenzodioxin targeted both CXCL12 and ICAM1 (Figure 7a). However, streptozocin is an antibiotic commonly used for diabetes modeling, dietary fats represent the collective oils from food sources, and tetrachlorodibenzodioxin is a highly toxic, persistent organic pollutant with no reported pharmacological application. Therefore, only dexamethasone was retained as the key compound for molecular docking. Docking analysis demonstrated that CXCL12 bound to dexamethasone with a binding energy of −6.2 kcal/mol, forming one hydrogen bond (bond length: 2.9 Å) (Figure 7b). ICAM1 bound to dexamethasone with a binding energy of −6.2 kcal/mol, forming two hydrogen bonds (bond lengths: 1.8 and 2.6 Å) (Figure 7c), suggesting a stable binding conformation with favorable affinity.

Additionally, the results of the MD simulation showed that the RMSD values of CXCL12 and ICAM1 were concentrated between 0.5 and 1, which suggested that their proteins were in a stable state from 1 to 50 ns and 5 to 100 ns (Figure 8a). Meanwhile, the RMSF values of CXCL12 and ICAM1 were concentrated in the ranges of 0.1–0.25 and 0.1–0.2, respectively; this indicated that the amino acid flexibility of these proteins was favorable and that their binding to the dexamethasone ligand was stable (Figure 8b). Furthermore, the energy values of CXCL12 and ICAM1 were low during the simulation and exhibited minimal fluctuations, which demonstrated that the binding of CXCL12 and ICAM1 to dexamethasone was relatively stable (Figure 8c).

Clinical samples were analyzed by RT-qPCR to validate the expression of ICAM1, CXCL12, and MMP3. The results confirmed that all three genes were significantly upregulated in PD samples (p < 0.001) (Figures 9a, 9b, and 9c), consistent with the expression profiles observed in transcriptome sequencing data. To further validate these findings at the protein level, IHC staining was performed on periodontal tissues from PD patients and healthy controls (n = 10 per group) using antibodies against ICAM1, CXCL12, and MMP3. As shown in Figure 9d, strong immunopositivity was observed for all three proteins in PD tissues compared with controls, with staining predominantly localized to epithelial and connective tissue regions of the gingiva.

Quantitative analysis of IOD further confirmed significantly elevated expression levels of ICAM1, CXCL12, and MMP3 in PD tissues relative to controls (p < 0.001) (Figures 9e, 9f, and 9g). These protein-level results were consistent with mRNA expression data from transcriptomic datasets and RT-qPCR validation, reinforcing the reliability of the findings. The upregulation of these SASP-related proteins in PD supports their potential role in mediating local inflammation and immune dysregulation.