Bioinformatics identification of forkhead box O3 and the proto-oncogene tyrosine-protein kinase receptor RET as aging-associated biomarkers linked to immune cell infiltration in Parkinson’s disease

Data collection and processing

The gene expression datasets GSE20163, GSE20164, and GSE8397 were obtained from the publicly available National Center for Biotechnology Information (NCBI) gene expression omnibus (GEO) database (https://www. ncbi.nlm.nih.gov/geo/). These three datasets consist of data derived from SN tissues, including 25 healthy controls (HCs) and 38 patients with PD, and were utilized as the test dataset. In addition, GSE7621 (9 HCs and 16 PD patients) and GSE49036 (8 HCs and 15 PD patients) served as the validation datasets. Furthermore, 307 aging-related genes were extracted from the human aging gene database (HAGR)^[21] (http://genomics.senescence.info/genes/).

Weighted gene co-expression network analysis (WGCNA)

WGCNA^[22] was employed to identify the relationships between gene modules and PD. By employing the R package “WGCNA,” weighted co-expression networks were constructed using genetic and clinical data from GSE20163, GSE20164, and GSE8397. The slope was set to approximately 1, and the scale-free soft threshold power R² was adjusted to approach 0.85 to convert the adjacency matrix into a topologically overlapping matrix. The adjacency matrix was then transformed into the topological overlap matrix (TOM) using an optimal power (β = 8). After the modules were identified, hierarchical clustering was applied to calculate the Eigen genes. Finally, Spearman’s correlation analysis was performed to evaluate the relationships among the different modules.

Functional enrichment analysis

Gene ontology (GO) and kyoto encyclopedia of genes and genomes (KEGG) enrichment analyses were conducted using the clusterProfiler package in R to identify the functions of key genes involved in PD.^[22] The analysis included biological processes, cellular components, and molecular functions, with statistical significance set at P < 0.05.

Protein–protein interaction (PPI) network

The online STRING platform (https://cn.string-db.org/) was utilized to construct the PPI network, and a combined interaction score threshold greater than 0.4 was applied. The network was visualized using Cytoscape.^[23]

Validation of diagnostic risk models

The predictive performance of aging-associated genes for identifying PD in both the training and validation cohorts was assessed using least absolute shrinkage and selection operator (LASSO) regression and ROC curve analyses. The diagnostic accuracy was evaluated based on the area under the curve (AUC), which quantified the ability to distinguish PD patients from HCs.

Gene set enrichment analysis (GSEA)

The potential pathways related to the identified diagnostic genes were examined through GSEA implemented using the R package GSEA. Significant enrichment was defined as an adjusted P < 0.05.

Immune cell infiltration abundance analysis

Immune cell infiltration in PD and normal tissues was quantified using the ImmuCellAI platform new version (http://bioinfo.life.hust.edu.cn/web/ImmuCellAI/), covering 24 immune cell populations, including 18 T-cell subsets and six other immune cell types (DCs, B cells, macrophages, monocytes, neutrophils, and NK cells).^[24] The correlations between the diagnostic genes and immune cell infiltration levels in PD patients were assessed through Spearman’s correlation analysis, and the results were visualized with an R-based heatmap.

Cell culture

Human neuroblastoma (SH-SY5Y) cells (SCSP-5014), purchased from the Shanghai Cell Bank of the Chinese Academy of Sciences, were cultured in DMEM (41401ES76, YEASEN, China) or DMEM/F-12 (41406ES76, YEASEN, China) supplemented with 10% fetal bovine serum (FSP500, EXCELL), 100 U/mL penicillin, and 100 μg/mL streptomycin in a humidified incubator at 37°C with 5% CO₂. Afterward, the cells were incubated with 500 µM 1-Methyl-4-phenylpyridinium (MPP⁺; M10041, AbMole, USA) to induce a PD cell model or with 200 µM H₂O₂ (Lircon, China) to induce oxidative damage in SH-SY5Y cells, as described in previous studies.^[25,26] The cell lines were authenticated by short-tandem repeat profiling, and tests for mycoplasma contamination were confirmed to be negative.

Extraction of RNA followed by quantitative real-time polymerase chain reaction (qRT-PCR)

Total ribonucleic acid (RNA) was isolated from SH-SY5Y cells using TRIzol reagent (15596026CN, Invitrogen, USA) according to the manufacturer’s guidelines. The concentration and purity of the extracted RNA were determined using a NanoDrop 2000 spectrophotometer (Thermo Scientific, USA). One microgram of RNA was then reverse transcribed into complementary DNA (cDNA) using RNA-to-cDNA Master Mix (4387406; Applied Biosystems™, Thermo Scientific, USA) according to the provided instructions.

qRT-PCR was performed on a Rotor-Gene Q instrument (Qiagen, Germany) utilizing SYBR Premix Ex Taq II (RR041A, TaKaRa, Japan) to quantify RET and forkhead box O3 (FOXO3) gene expression. The thermal cycling conditions were as follows: initial denaturation at 95°C for 15 min, followed by 45 amplification cycles of 95°C for 5 s and 60°C for 30 s. Rotor-Gene Real-Time Analysis Software 6.0 was used for data analysis. Gene expression levels were normalized to those of the internal control GAPDH, and relative expression was calculated using the 2^−ΔΔCt method. The following primers were used:

• RET-forward 5’TCCTTTCCCTTACCCACCTTCAG3’,

• RET-reverse 5’AGACCACAGCACCACAGACC3’,

• FOXO3-forward 5’GCAGACCATCCAAGAGAACA AGC3’,

• FOXO3-reverse 5’TGGCTAAGTGAGTCCGAAGTGAG3’,

• GAPDH-forward 5’TCAAGGCTGAGAACGGGAAG3’,

• GAPDH-reverse 5’CGCCCCACTTGATTTTGGAG3 ’.

In vivo experimental study

A PD animal model was established by intraperitoneally injecting 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine (MPTP, 30 mg/kg; M0896, Sigma-Aldrich, USA) into C57BL/6 mice (age: 6-8 weeks, weight: 25 ± 2 g) for five consecutive days. Twenty-month-old aging mice on a C57BL/6J background were obtained from Phenotek Biotechnology Co., Ltd. (Shanghai, China). All mice were maintained in a specific pathogen-free facility under standardized conditions, including a constant temperature of 23 ± 2°C, a 12-h light/dark cycle, and a relative humidity of 60% ± 5%. Throughout the experiment, the animals had unrestricted access to standard chow and water. The mice were euthanized at the completion of the experiment using an intraperitoneal injection of pentobarbital sodium (110 mg/kg; P3761, Sigma, Merck, USA). All procedures involving animals were performed in compliance with the ARRIVE guidelines and Guidelines for the Care and Use of Laboratory Animals in China and received approval from the Experimental Animal Center of Tongji University School of Medicine (Approval No. 2023-MHFY-30JZS).

Immunohistochemistry

Following deep anesthesia, transcardial perfusion was performed using phosphate-buffered saline (PBS), followed by 4% paraformaldehyde (PFA) in PBS. Afterward, the brains were fixed with 4% PFA for 24 h and sliced into 30 μm thick sections. The sections were incubated with an anti-RET antibody (ab134100, Abcam, UK), an anti-FOXO3 antibody (AFRM9032, AiFang Biological, China), or a tyrosine hydroxylase antibody (ab6211, Abcam, UK). The experiments followed the standard protocol, and the images were captured with an optical microscope (Axio Observer 7, Zeiss, Germany).

Mitochondrial ROS measurement

Mitochondrial ROS levels were detected using MitoSOX Red (S0061, Beyotime, China). The cells were incubated with a working solution containing MitoTracker (C1048, Beyotime, China; diluted 1:5000) and MitoSOX (1:1000 dilution in PBS) for 30 min at 37°C, washed with PBS, and then subjected to confocal microscopy (Olympus FV1200) for visualization.

Statistical analysis

Statistical analyses were conducted using R software (version 4.0.3) and GraphPad Prism (version 9.0). Spearman’s correlation analysis was performed to assess the relationships between variables. The data are presented as the means ± standard deviations (SDs), for comparisons between two groups, unpaired Student’s t-tests were employed, whereas one-way analysis of variance (ANOVA) was used to evaluate differences among multiple groups. When significant differences were detected by ANOVA, Tukey’s post hoc test was performed for pairwise comparisons. P <0.05 was considered to indicate statistical significance.

Research design flow chart

This flowchart outlines the process used to analyze gene expression [Figure 1]. The training dataset included 25 HCs and 38 PD patients from GSE20163, GSE20164, and GSE8397. The process began with aging-related differentially expressed genes (DEGs) identified through WGCNA. These genes were subjected to GO and KEGG enrichment analyses and PPI network construction. LASSO regression was subsequently applied, followed by a ROC curve analysis to determine the hub aging-related DEGs. These genes were subsequently validated using the GSE7621 (9 HCs and 16 PD patients) and GSE49036 (8 HCs and 15 PD patients) datasets. In parallel, an immune cell infiltration analysis was conducted using the ImmuCellAI tool, followed by Spearman’s correlation analysis. GSEA was performed to further explore the biological significance of the hub aging-related DEGs.

Close

Figure 1: Bioinformatics workflow for identifying biomarkers related to aging and immune cell infiltration in PD patients. HC: Healthy control, PD: Parkinson’s disease, WGCNA: Weighted gene co-expression network analysis, DEGs: Differentially expressed genes, GO: Gene ontology, KEGG: Kyoto encyclopedia of genes and genomes, LASSO: Least absolute shrinkage and selection operator, PPI: Protein‒protein interaction, ROC: Receiver operating characteristic, GSEA: Gene set enrichment analysis.

Export to PPT

Identification of aging-related genes in PD patients

First, the WGCNA network was constructed using the GSE20163, GSE20164, and GSE8397 datasets. A scale-free network was then constructed with a β value of 8 (R² = 0.85) [Figure 2a and b]. Twenty-six co-expression modules were subsequently identified [Figure 2c]. The correlation coefficients of each module with the PD group were calculated to further investigate the relationships between the modules and the phenotype. The results revealed that the turquoise module (r = 0.35, P = 1.4e−14) and the light cyan module (r = 0.77, P = 1.2e−40) exhibited significant positive correlations with the PD group [Figure 2d and e]. A total of 656 genes from the turquoise and light cyan modules were selected for further analysis.

Close

Figure 2: Screen of Parkinson’s disease aging-related genes. (a) Network connectivity under different soft thresholds in Weighted gene co-expression network analysis (WGCNA). (b) Scaleless network display of the optimal soft threshold in WGCNA. (c) Correlation analysis results between the gene clustering module and PD. (d) Co-expression correlations among gene members in the turquoise gene module. (e) Co-expression correlation of gene members within the light cyan gene module.

Export to PPT

Enrichment analysis of aging-related genes and PPI network analysis

We initially downloaded 307 aging-related genes from the aging database. We subsequently intersected the 656 module genes with the 307 aging-related genes to identify 18 aging-related genes significantly associated with the aging group [Figure 3a]. GO functional classification and KEGG enrichment analyses were performed to elucidate the molecular roles of these parkinson aging-related genes (PARGs). The cellular component analysis revealed that PARGs were predominantly localized in the cytosol, protein-containing complexes, and membrane protein complexes [Figure 3b]. The biological processes in which these genes were involved included the response to organic substances, cellular response to organic substances, intracellular signal transduction, and response to oxygen-containing compounds [Figure 3c]. The molecular function analysis indicated that the aging-related genes were linked to enzyme binding, signaling receptor binding, protein kinase binding, and kinase binding [Figure 3d]. Through the KEGG enrichment analysis, aging-associated genes were shown to participate in key pathways, including PI3K–Akt signaling, neurotrophin signaling, AMPK signaling, and pathways regulating longevity [Figure 3e]. A PPI network was constructed for the 18 aging-related genes, which included protein phosphatase, Mg²⁺/Mn²⁺ dependent 1D (PPM1D), breast cancer susceptibility gene 1(BRCA1), small ubiquitin like modifier 1 (SUMO1), protein-L-isoaspartate (D-Aspartate) O-methyltransferase (PCMT1), eukaryotic translation elongation factor 1 epsilon 1 (EEF1E1), angiotensin II receptor type 1 (AGTR1), RET, FOXO3, succinate dehydrogenase complex subunit C (SDHC), CAMP responsive element binding protein 1 (CREB1), insulin receptor (INSR), insulin receptor substrate 1 (IRS1), suppressor of cytokine signaling 2 (SOCS2), growth hormone receptor (GHR), apoptosis inducing factor mitochondria associated 1 (AIFM1), mitogen-activated protein kinase Kinase Kinase 5 (MAP3K5), cell division cycle 42 (CDC42), and phospholipase C gamma 2 (PLCG2) [Figure 3f]. These findings suggest a significant imbalance in aging-related genes and pathways in PD patients.

Close

Figure 3: Enrichment analysis and construction of a PPI network of aging-related genes. (a) Eighteen genes were obtained by crossing the differential genes and aging-related genes obtained from the WGCNA. (b-d) GO functional enrichment analysis (cellular components, biological processes, and molecular functions). (e) KEGG pathway enrichment analysis. (f) PPI network of PARGs. PPI: Protein‒protein interaction, WGCNA: Weighted gene co-expression network analysis, PARGs: Parkinson aging-related genes.

Export to PPT

Identification of optimal biomarkers in PD and analysis of the diagnostic risk model

The results are illustrated in LASSO coefficient trajectory plots [Figure 4a] and LASSO regression model diagrams [Figure 4b]. The results indicated that the LASSO regression model included seven genes: BRCA1, CREB1, FOXO3, INSR, RET, SDHC, and SUMO1. In addition, ROC curves were plotted [Figure 4c-e] for both the LASSO regression model and PD samples from the GSE7621 and GSE49036 datasets. The findings demonstrated that the PD diagnostic risk model exhibited high accuracy in both the test set (AUC = 96.3158) and the two validation sets, GSE7621 (AUC = 79.8611) and GSE49036 (AUC = 77.5000).

Close

Figure 4: LASSO regression and analysis of the diagnostic model. (a) Variable trajectory diagram. (b) Diagram of the diagnostic risk model from the LASSO regression model. (c) The ROC curve of the diagnostic model. (d) The ROC curve of the samples in GSE7621. (e) The ROC curve of the samples from GSE49036. LASSO: Least absolute shrinkage and selection operator, ROC: Receiver operating characteristic.

Export to PPT

Verification of the expression levels of optimal biomarkers

Differences in the expression of optimal biomarkers between HCs and PD patients were further examined in two independent validation datasets (GSE7621 and GSE49036) to assess the roles of key aging-associated genes [Figure 5a-n]. The expression level of FOXO3 was increased [Figure 5c and j], whereas that of RET [Figure 5e and l] was significantly decreased in the PD group compared with the HC group across both datasets (GSE7621 and GSE49036). GSEA was employed to elucidate the relationships between gene expression and the biological processes, cellular components, and molecular functions affected in PD patients and to explore the overall impact of gene expression in PD samples [Figure 5o and p]. The findings indicated that the genes were significantly enriched in biological functions and signaling pathways associated with oxidative physiology, the tricarboxylic acid cycle (TCA) amino sugar, and nucleoside sugar metabolism and PD, but tended to be associated with cytokine and cytokine receptor interactions, the MAPK signaling pathway, linoleic acid metabolism, and the JAKSTAT signaling pathway [Figure 5o and p]. These results suggest that these genes play crucial roles in the development of PD. In addition, we observed that all the genes were associated with immune-related pathways.

Close

Figure 5: Verification of the expression levels of aging-related genes and GSEA of aging-related genes. (a-g) In GSE7621, BRCA1 (a), CREB1 (b), FOXO3 (c), INSR (d), RET (e), SDHC (f), and SUMO1 (g) expression differed between HCs and PD patients. (h-n) In GSE49036, BRCA1 (h), CREB1 (i), FOXO3 (j), INSR (k), RET (l), SDHC (m), and SUMO1 (n) expression differed between HCs and PD patients. (o) GSEA showed that the module genes significantly affected oxidative physiology (P = 0, FDR = 0.0016), the TCA cycle (P = 0.002, FDR = 0.0087), amino sugar and nuclear sugar metabolism (P = 0.002, FDR = 0.014), and Parkinson’s disease (P = 0.0022, FDR = 0.0032). (p) GSEA showed that the module genes significantly affected cytokine and cytokine receptor interaction (P = 0.0078, FDR = 1), the MAPK signaling pathway (P = 0.0099, FDR = 0.67), linoleic acid metabolism (P = 0.012, FDR = 1), and the JAK-STAT signaling pathway (P = 0.012; FDR = 0.77). The screening criterion for the GSEA was an adjusted P < 0.05. FDR: False discovery rate, FOXO3: Forkhead box O3, HC: Healthy control, PD: Parkinson’s disease, GSEA: Gene set enrichment analysis.

Export to PPT

Immune cell infiltration analysis

An immune cell infiltration analysis was performed using the expression matrix of PD and HC samples from the test datasets (GSE20163, GSE20164, and GSE8397). The abundance of 24 immune cells in both the PD and HC groups was calculated using the single-sample Gene Set Enrichment Analysis (ssGSEA) algorithm. A violin plot comparing the groups [Figure 6a-c] revealed significant differences in the levels of five types of immune cells, namely, monocytes, macrophages, natural killer T cells, naive CD4⁺ cells, and Tr1 cells (P ≤0.05). A correlation heatmap was subsequently generated to illustrate the correlation results of the immune cell infiltration abundances in the PD group [Figure 6d]. Within the samples from the PD group, the numbers of NKT cells were negatively correlated with monocyte (correlation coefficient = -0.38), naive CD4⁺ T-cell (correlation coefficient = -0.35), and Tr1 (correlation coefficient = -0.41) counts. In addition, naive CD4⁺ T cells were positively correlated with monocytes (correlation coefficient = 0.47) and Tr1 cells (correlation coefficient = 0.60).

Close

Figure 6: Analysis of immune cell infiltration using the ImmuCellAI. (a) Comparative violin plot of the immune cell infiltration scores of PD patients and HCs. (b and c). Violin plots comparing the abundance of immune cell infiltration in HCs and PD patients. (d) The results of the correlation analysis of monocytes, macrophages, NKT cells, naive CD4⁺ cells, and Tr1 cells in the PD group. Significance levels were denoted as follows: ^✶P < 0.05, ^✶✶P < 0.002, ^✶✶✶P < 0.0002; ^✶✶✶✶P < 0.0001. HCs: Healthy controls, PD: Parkinson’s disease.

Export to PPT

Correlations between aging-related genes and immune cells

We observed significant differences in the expression levels of two aging-related genes, RET and FOXO3, between the PD and HC groups in both the GSE7621 and GSE49036 datasets. Therefore, we focused on these two genes for further analysis. An analysis of immune cell infiltration revealed significant differences in the abundance of the five immune cell types between the two groups, with all the differences reaching statistical significance. As shown in Figure 7, the results of Spearman’s correlation analysis indicated both positive and negative relationships. RET was negatively correlated with naive CD4⁺ T cells (R = -0.344, P = 0.006), macrophages (R = -0.290, P = 0.021), monocytes (R = -0.333, P = 0.008), and Tr1 cells (R = -0.458, P <0.001), whereas it was significantly positively correlated with NKT cells (R = 0.281, P = 0.026) [Figure 7a-e]. FOXO3 was not significantly correlated with naive CD4⁺ T cells, macrophages, monocytes, or NKT cells [Figure 7f-i] but was significantly positively correlated with Tr1 cells (R = 0.311, P = 0.013) [Figure 7j].

Close

Figure 7: Correlations between aging-related genes and immune cells. (a-e) Diagrams showing the correlations of the RET expression level with the abundance of infiltrating naive CD4⁺ T cells (a), macrophages (b), monocytes (c), NKT cells (d), and Tr1 cells (e). (f-j). Diagrams showing the correlations of the Forkhead box O3 expression level with the abundance of infiltrating naive CD4⁺ T cells (f), macrophages (g), monocytes (h), NKT cells (i), and Tr1 cells (j).

Export to PPT

Verification of the expression of aging-related genes in in vivo and in vitro models

RNA expression was analyzed in cellular models of PD and aging, established through MPP⁺ or H₂O₂ treatment, to assess the feasibility of using aging-associated biomarker genes in PD diagnosis. MPP⁺ is taken up by neurons and accumulates in mitochondria, where it inhibits complex I of the electron transport chain. This inhibition leads to the production of reactive oxygen species (ROS), ultimately resulting in neuronal damage.^[26] Costaining with MitoTracker and MitoSOX demonstrated robust mitochondrial ROS production induced by MPP⁺ in SHSY5Y cells, confirming the establishment of an in vitro PD model [Supplementary Figure 1a and b]. The qRT-PCR results showed that treatment with MPP⁺ or H₂O₂, which was used to induce PD or aging, led to the significant upregulation of FOXO3 expression and downregulation of RET expression in the cells [Figure 8a and b]. MPTP was used to establish an in vivo PD model, which caused a significant loss of tyrosine hydroxylase (TH)-positive neurons in the SN, a hallmark pathological feature of PD [Supplementary Figure 1c and d]. As shown in [Figure 8c and d], moderate FOXO3 immunoreactivity was detected in the SN of control mice. On MPTP treatment, which was used to induce the PD model, a marked increase in the FOXO3 staining intensity was observed in aged mice, indicating elevated levels of FOXO3 in both pathological and physiological aging states. In contrast, RET immunoreactivity was prominently expressed in the control group but dramatically decreased following MPTP treatment. Compared with control mice, aged mice also exhibited reduced RET expression, although the decrease was less pronounced in aged mice than in MPTP-treated mice [Figure 8e and f]. These findings suggest that these two aging-related genes may serve as reliable diagnostic biomarkers for PD.

Close

Figure 8: Forkhead box O3 (FOXO3) and RET expression in in vivo and in vitro models. (a and b) Quantitative real-time polymerase chain reaction analysis showing the relative mRNA expression levels of FOXO3 (a) and RET (b) in SH-SY5Y cells treated with MPP⁺ (500 µM) or H2O2 (200 µM) for 24 h compared with those in normal control (NC) cells. n = 3 independent experiments. (c) Representative images of immunohistochemical (IHC) staining for FOXO3 in the SN of mice from the control, MPTP-treated, and aging groups. Scale bar = 20 µm. (d) Quantification of the FOXO3 staining intensity in (c), normalized to that of the NC. n = 3 mice per group. (e) Representative images of IHC staining of RET in the SN across the same three groups. Scale bar = 20 µm. (f) Quantification of the RET staining intensity in (e), normalized to that of the NC. n = 3 mice per group. The data are presented as the means ± SDs. Statistical analysis was performed using one-way ANOVA followed by Tukey’s post hoc test. ^✶P < 0.05, ^✶✶✶P < 0.0002, ^✶✶✶✶P < 0.0001, ns: Not significant. NC: Normal control, MPP+: 1-Methyl-4-phenylpyridinium, MPTP: 1-methyl-4-phenyl-1,2,3,6-tetrahydropyridine, H₂O₂: hydrogen peroxide, FOXO3: Forkhead box O3, RET: Proto-oncogene tyrosine-protein kinase receptor RET.

Export to PPT