Identification and Experimental Validation of Parkinson’s Disease with Major Depressive Disorder Common Genes

Parkinson’s disease (PD) is the second most common neurodegenerative disease that affects about 10 million people worldwide. Non-motor and motor symptoms usually accompany PD. Major depressive disorder (MDD) is one of the non-motor manifestations of PD it remains unrecognized and undertreated effectively. MDD in PD has complicated pathophysiologies and remains unclear. The study aimed to explore the candidate genes and molecular mechanisms of PD with MDD. PD (GSE6613) and MDD (GSE98793) gene expression profiles were downloaded from Gene Expression Omnibus (GEO). Above all, the data of the two datasets were standardized separately, and differentially expressed genes (DEGs) were obtained by using the Limma package of R. Take the intersection of the two differential genes and remove the genes with inconsistent expression trends. Subsequently, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analyses were investigated to explore the function of the common DEGs. Additionally, the construction of the protein–protein interaction (PPI) network was to search the hub genes, and then the least absolute shrinkage and selection operator (LASSO) regression was used to further identify the key genes. GSE99039 for PD and GSE201332 for MDD were performed to validate the hub genes by the violin plot and receiver operating characteristic (ROC) curve. Last but not least, immune cell dysregulation in PD was investigated by immune cell infiltration. As a result, a total of 45 common genes with the same trend. Functional analysis revealed that they were enriched in neutrophil degranulation, secretory granule membrane, and leukocyte activation. LASSO was performed on 8 candidate hub genes after CytoHubba filtered 14 node genes. Finally, AQP9, SPI1, and RPH3A were validated by GSE99039 and GSE201332. Additionally, the three genes were also detected by the qPCR in vivo model and all increased compared to the control. The co-occurrence of PD and MDD can be attributed to AQP9, SPI1, and RPH3A genes. Neutrophils and monocyte infiltration play important roles in the development of PD and MDD. Novel insights may be gained from the findings for the study of mechanisms.


Introduction
Parkinson's disease is the second most common neurodegenerative disorder with an increasing incidence rate worldwide [1]. PD is characterized by age-related loss of dopaminecontaining neurons in the substantia nigra pars compacta (SNpc) region and the formation of Lewy bodies which consisted of α-synuclein [2]. These pathological characteristics manifested with motor symptoms of bradykinesia, rigidity, postural instability, and static tremor, as well as non-motor symptoms, such as sleep disturbances, depression, hyposmia, and constipation [3]. Studies have reported that non-motor symptoms are earlier and more serious than motor symptoms [4,5]. Therefore, it is important to pay early attention to the non-motor symptoms.
Depression (major depressive disorder, MDD) is one of the non-motor symptoms and is considered the most common neuropsychiatric disorder with PD. An increase in death risk is associated with MDD due to poor cognitive performance, poor quality of life, and a worse functional status [6]. Studies have demonstrated that patients with PD are more susceptible to MDD compared with healthy controls [7]. Cross-sectional studies have found that 2.7 to 90% of PD patients have MDD [8]. It is crucial to diagnose PD patients early since MDD significantly impacts their quality of life. However, it is challenging to diagnose timely and accurate MDD in patients with PD, because many features of depression, such as weight loss, insomnia, psychomotor retardation, and fatigue, overlap with the primary symptoms of PD or with the side effects of medication [8].
There are several theories about how MDD occurs in PD, although the exact mechanism remains unclear [9]. It is widely accepted that MDD is caused by monoamine oxidase. PD patients with MDD postmortem had substantial destruction in the noradrenaline-producing nucleus coeruleus of the brain. Additionally, morphological changes in the nucleus raphe were found in PD patients with MDD [10]. It appears that MDD in PD is caused by the degeneration of monoaminergic neurotransmitter systems and dysfunction of the fronto-cortical system. Besides, the dysfunction of the orbitofrontal cortex is induced by the degeneration of dopaminergic mesocortical and mesolimbic neurons, which affects serotonergic cell bodies in the dorsal raphe nuclei secondarily [11]. Neuroinflammation is another aspect that affects the production of serotonin which is linked to MDD and neurodegeneration [12]. Above all, the distinctive mechanisms that explain the cooccurrence of PD and MDD remain unclear.
The current investigation examined the transcriptome expression data of whole blood (GSE98793 and GSE6613) obtained from patients diagnosed with MDD and PD in the GEO database to investigate their shared molecular mechanisms. Additionally, gene modules were analyzed by constructing protein-protein interaction (PPI) nodes and utilizing the search tool for the retrieval of interacting genes/proteins (STRING) database and Cytoscape software to identify hub genes. LASSO regression was employed to further identify the crucial genes. Finally, the hub genes were validated using GSE99039 (PD) and GSE201332 (MDD). Based on an analysis of independent datasets (GSE99039), we verified the expression level of hub genes and looked at the correlation of PD clinical characteristics. Besides, we also constructed models to further validate the hub genes by qPCR and WB. The results of this study might provide new insights into PD and MDD pathogenesis mechanisms.

Data Collection and Preprocessing
GSE6613 [13,14] and GSE99039 [15]) of PD datasets as well as (GSE98793) [16] and (GSE201332) [17] of MDD datasets were downloaded from the GEO (https:// www. ncbi. nlm. nih. gov/ geo/) database [18]. Detailed information on the 4 datasets is shown in Table 1, and a flow chart of the study design is shown in Fig. 1. The R software 4.0.5 was used to normalize, log2 transform, and convert gene names for each dataset.

Identification of Differentially Expressed Genes
Microarray data were downloaded from GSE6613 and GSE98793, and box plots and PCA analysis were used to represent the normalized expression matrix. Multiple probes identifying the same gene were averaged to determine its expression level. |log2 Fold change (FC)|> 0 (PD filtration) or |log2 Fold change (FC)|> 0.585 (MDD filtration) and p value < 0.05 were set as the criteria for identifying DEGs using the "Limma" package of R software [19]. A variety of heat maps, volcano maps, and box line plots were created using the "heatmap" and "ggplot2" packages of R software (version 3.6.3) [20].

Protein-Protein Interaction Network Construction
Based on the identified common DEGs with the same expression trend, STRING (https:// string-db. org/) is a database, and used to construct the protein-protein interaction (PPI) network with combined scores greater than 0.4 [23]. The PPI network was then constructed and presented on the Cytoscape platform [24]. The significant modules and core genes were identified using CytoHubba, a plugin in Cytoscape. Four different algorithms (Maximum Neighborhood Component (MNC), Maximal Clique Centrality (MCC), DEGREE, and Density of Maximum Neighborhood Component (DNMC)) were used to identify hub genes. Lastly, a Venn diagram was used to demonstrate that the common genes obtained by the four algorithms were reliable hub genes.

Identification of Common Hub Genes Using LASSO Logistic Regression
LASSO analysis is a regression method that improves prediction accuracy by selecting a variable from high dimensional data that has a strong predictive value and low correlation [25]. The hub genes were further identified by LASSO analysis. LASSO logistic regression was built using hub gene expression levels and clinical traits. To differentiate patients from healthy controls, ROC curves were used using the "pROC" R package [26]. Moreover, the hub genes were also represented by box line plots.

Validation of Hub Gene Expression
GSE99039 for PD [15] and GSE201332 for MDD [17] were used to validate the expression levels of the identified hub genes. GSE99039 contains 205 PD and 233 control samples. GSE201332 contains 20 MDD and 20 control samples. The expression difference of hub genes was also represented by a box line plot and ROC curve. Wilcoxon test was used to compare the two datasets. p value < 0.05 was considered significant.

Animals
Male C57BL/6 mice weighing 23-27 g at the age of 9-11 weeks were purchased from Pengyue experimental animal Ltd. (Jinan, Shandong, People's Republic of China) and acclimated for at least 1 week. All mice were housed under controlled lighting conditions in a 12/12-h light/dark cycle and a comfortable temperature of 22-26 °C with water available ad libitum but the food required our experiments. Additionally, all animal care and experimental procedures were approved by the Ethics Committee. These animals' conditions are consistent with our previous studies [27].

Chronic Unpredictable Mild Stress Model
The model of depression was processed with chronic unpredictable mild stress (CUMS), the protocol was described as the previous study [29]. The animals in various groups were unpredictably subjected to mild stress for 4 weeks except for those animals in the vehicle + no stress group. Individual stressors and the length of time they were applied each day were as follows: (i) food deprivation; (ii) water deprivation; (iii) restraint; (iv) restraint at 4•C; (v) flashing light for 120-210 min; and (vi) isolation. Stressor stimuli were applied at different times every day, to minimize their predictability.

Immunohistochemistry
Frozen sections of the SNpc were conducted as the protocol of the immunohistochemistry kit as our previous study described [27]. The number of TH-positive neurons represented dopaminergic degeneration.

Rota-Rod Training
The mice were trained on the Rota-rod apparatus (Anhui Zhenghua Biological Apparatus Company, Anhui, China) with 15 rpm for 20 min/day for 1 week. After pretraining, animals were assigned to two groups (n = 6 per group) treat with or not MPTP.

Pole Test
An experimental apparatus consisting of a wooden rod with a diameter of 1 cm and a height of 50 cm, wrapped in non-adhesive gauze and topped with a wooden ball with a diameter of 2 cm, was vertically positioned within a cage containing bedding material. The experimental procedure involved placing a mouse in a head-down position on the ball and recording the time taken to climb from the ball to the bottom of the cage. Each mouse was subjected to three trials with an interval of one hour, and the mean value was subsequently calculated for further analysis.

Forced Swimming Test
Forced swimming test (FST) was conducted followed by the previous study [30]. In this study, the FST was carried out to evaluate depressive-like behavior associated with stress or MPTP effects. The FST apparatus consisted of a glass cylinder having a 25 cm diameter and a 35 cm height. Water was maintained at a temperature between 22 and 26 °C. The various groups were subject to FST for 6 min to evaluate the immobility, climbing, and swimming times.

RNA Isolation and Quantitative Polymerase Chain Reaction
Total RNA was extracted from the peripheral blood of multiple groups according to the manufacturer's protocol (Trizol method and miRNeasy Mini kit (Qiagen)). To reverse transcribe total RNA, the HiScript Q RT Super-Mix for quantitative polymerase chain reaction (qPCR) Kit (Vazyme, R123-01) was used. In this experiment, QuantStudio 6 (Applied Biosystems) was used to perform quantitative PCR to measure SYBR Green levels in realtime PCR using the manufacturer's recommended cycling conditions. In the analysis, duplicates of samples with Cq values over 35 were removed. The primers were provided in Table 1. The 2 − ∆∆Ct method was applied to calculate the relative expression level of mRNA.

Western Blot Analysis
Equal quantities of the protein (30 μg) were subjected to SDS-PAGE (4-20% gel) and transferred to a PVDF membrane. After blocking in 5% non-fat dry milk in TBS containing Tween 20 (TBST) for 1 h at room temperature, Subsequently, the membranes were incubated with specific primary antibodies at 4 °C overnight. The following day, the membranes were washed three times with TBST and subsequently incubated with the corresponding secondary antibody for 1 h at room temperature. Then, the membranes were washed again with TBST 3 times, 10 min each time. The epitope was visualized by an enhanced chemiluminescence system (New Cell and Molecular Biotech, China). The gray value of the target band was analyzed by Image J software.

Gene Set Enrichment Analysis
Based on gene expression levels for population phenotypes, GSEA (http:// softw are. broad insti tute. org/ gsea/ index. jsp) was used to assess pathway and molecular mechanisms relationships between the two groups [31]. Enriched gene sets with nominal p values of < 0.05, |normalized enrichment scores (NES) |> 1, and false positive rate (FDR) q values of < 0.25 were considered significant.

Immune Infiltration Analysis
CIBERSORT, a computational method for estimating the proportion of diverse immune cells based on gene expression profiles, was used to assess immune cell proportion in patients and healthy control [32]. Infiltration analysis using the "Cibersort" R package was performed. Bar plots were used to visualize the proportion of different types of immune cells. Box plots were used to compare the proportions of different types of immune cells in two datasets. The "corrplot" R package was used to create a heatmap illustrating the correlation between 22 types of infiltrating immune cells [33].

Statistical Analysis
All statistical analyses were performed using R software. P value < 0.05 was considered statistically significant.

Identification of Common DEGs Between PD and MDD
The expression matrices of the two datasets GSE6613 and GSE98793 were normalized, and box plots showed straight lines as distribution trends ( Fig. 2A-D). Data repeatability between the two datasets was assessed using PCA in this study, and the results showed good repeatability between the two datasets ( Fig. 2E-F). For the GSE6613 dataset, according to the Limma R method, 797 DEGs were identified in the combined PD dataset, of which 420 were upregulated and 377 were downregulated, and a total of 2145 DEGs were identified as DEGs in MDD including 1384 were upregulated and 761 were downregulated. As shown in Fig. 3A, the volcano plot of PD DEGs was shown. Furthermore, the top 50 upregulated concurrently with the top 50 downregulated DEGs were shown using heatmaps (Fig. 3B, C). Regarding the MDD dataset, the volcano plot and heatmaps were shown in Fig. 3D-F. Besides, the intersection of the Venn diagram yielded 95 common DEGs (Fig. 4A). However, only 45 common DEGs with the same trend. The heatmaps of 45 common DEGs in the PD and MDD datasets were shown in Fig. 4B, C.

Enrichment Analysis of Common DEGs
To investigate the biological functions of the common DEGs, we performed KEGG pathways and GO terms enrichment analyses. Based on the KEGG enrichment analysis, common genes were primarily enriched in the "Metabolic signaling pathway," "Osteoclast differentiation," and "Steroid biosynthesis" (Additional Fig. 10A-B). Furthermore, GO analysis showed that genes were enriched in "neutrophil degranulation," "neutrophil activation involved in immune response," and "B cell activation involved in immune response" (BP); "cytoplasmic vesicle lumen," "secretory granule membrane," and "ficolin-1 rich granule" (CC); "immune receptor activity," guanyl-nucleotide exchange factor activity," and "protein phosphatase binding" (MF) (Additional Fig. 1C-F). As a result of these findings, the progression of both diseases appears to be mediated by inflammatory pathways as well as metabolism pathways.

Identification and External Validation of Candidate Hub Genes
Identifying hub genes that may be involved in the cooccurrence of PD and MDD. First of all, STRING was used to analyze PPI networks of communal DEGs to clarify interactions between them (Additional Fig. 2A). Modular gene enrichment analysis and significant gene modules (Additional Fig. 2B). MCC, MNC, Degree, and DNMC were used to predict and explore the top 15 hub genes in the PPI network using CytoHubba. 14 candidate hub genes were identified from the intersection of the 15 genes from the four algorithms: S100A11, COASY, SNCA, FECH, UGP2, CXCR1, AQP9, CSF3R, FDFT1, CFP, LILRB2, RPH3A, CYTH4, and SPI1 (Additional Fig. 2C). Moreover, GO and KEGG enrichment analysis was performed on the 14 hub genes (Additional Fig. 2D-E).
Subsequentially, nomogram construction and diagnostic value evaluation were performed using LASSO logistic regression. As shown in Fig. 5A, B, eight potential biomarkers were identified using the LASSO regression algorithm including CYTH4, AQP9, SPI1, CFP, RPH3A, FECH, UGP2, and FDFT1. Multivariate Cox regression analysis revealed that this prognostic model was an independent prognostic parameter of PD (Fig. 5C). Additionally, the heatmap showed the expression of the candidate genes according to the risk score (Fig. 5D). Through the LASSO analysis, the candidate genes were further validated. Fig. 3 The expression and identification of differentially expressed genes in the two datasets. A Volcano map of differentially DEGs between in GSE6613. The threshold was set to |log2FC (fold change) |> 0, and p-value < 0.05. B, C The top 50 upregulated concurrently with the top 50 downregulated DEGs of the PD dataset (D) volcano map of differential DEGs in GSE98793. The threshold was set to |log2FC (fold change) |> 0.585, and p-value < 0.05. E-F The top 50 upregulated concurrently with the top 50 downregulated DEGs of the MDD dataset ◂ At last, the candidate genes were further validated by external datasets. The expression of the screened eight genes was tested with GSE99039 of PD to improve the dependability of hub genes. As shown in Fig. 6A, AQP9, CFP, RPH3A, CYTH4, and SPI1 were consistent with the above study. In these three datasets, ROC curves were drawn using the R package with the expression of the four hub genes to assess the diagnostic accuracy. Five candidate genes possess a high diagnostic value (Fig. 6B-D). Additionally, we further applied the external validation of the MDD dataset with GSE201332. The violin plot showed that the expression of AQP9, RPH3A, and SPI1 was differentially expressed in the MDD dataset (Fig. 6E). Meanwhile, the ROC curve of three genes was established to have a significant value for diagnosing MDD (Fig. 6F). Above all, these results suggested that AQP9, RPH3A, and SPI1 could be promising markers for diagnosing PD and MDD.

In Vivo Model and FurtherValidated theHub Genes
Similar to some previous studies [27,34], MPTP-induced neurotoxicity of dopaminergic neurons exhibited as TH-positive neurons. As shown in Fig. 7A and Additional Fig. 3A, the MPTP group significantly showed the loss of TH-positive neuron populations in the SNpc and the depletion of dopaminergic nerve fibers in the striatum (STR) compared to the sham group. Additionally, the MPTP associated with the stress group also showed the same phenomenon. Besides, Rota-rod training and pole test all indicated the movement disorder of the PD model (Additional Fig. 3B, C). In other words, it also indicated the success of the PD model in vivo. As shown Fig. 7B-D is representing data on immobility, climbing, and swimming times in the forced swimming test. In the stress groups, the time of climbing and swimming was significantly decreased while the immobility time was enhanced as well as interaction for stress and MPTP.
The transcriptional changes of overlapped hub genes AQP9, RPH3A, and SPI1 were detected in the peripheral blood from the MPTP associated with stress groups and control groups by quantitative qPCR. The results indicated that the expression levels of AQP9, RPH3A, and SPI were both increased in the MPTP associated with stress groups in comparison with those in controls ( Fig. 8A-C), which was in line with the bioinformatics analysis. Besides, we also measured the protein expression in the animal model of AQP9, RPH3A, and SPI1. Interestingly, the expression of AQP9 and SPI1 were both increased in the MPTP associated with stress groups in comparison with those in the sham while the expression of RPH3A was decreased compared with the sham (Additional Fig. 4A-D).

GSEA Results of Hub Genes
Based on GSEA analysis of GSE6613 and GSE98793, differentially regulated pathways between the high and low expression groups were identified to determine the potential functions of hub genes in PD and MDD. GSEA analysis was conducted for the gene sets "PANTOTHENATE_AND_ COA_BIOSYNTHESIS," "GRAFT_VERSUS_HOST_DIS-EASE," "PROTEASOME," "CHRONIC_MYELOID_LEU-KEMIA," "NEUROTROPHIN_SIGNALING_PATHWAY" (Fig. 10).

Construction of Putative Hub Genes Protein-Protein Interaction Network and TF-Gene Interaction Network
By employing the tool "GeneMANIA," we constructed a putative PPI network of 20 genes involving hub genes APQ9, RPH3A, and SPI1 (Fig. 11A). Inner circle contains hub genes, while the outer circle contains predicted genes. We visualized the interaction network between candidate genes and TF genes using NetworkAnalyst (http:// www. netwo rkAna lyx. ca/). Human-mouse TF target interaction data are currently managed manually in the TRRUST database (https:// www. grnpe dia. org/ trrust/). Submitting target genes can be used to search for key TF (Fig. 11B).

Immune Cell Infiltration Analysis
Microenvironmental factors include immune cells, extracellular matrix, inflammatory factors, and growth factors that influence clinical therapeutic sensitivity and disease diagnosis. The CIBERSORT algorithm was used to estimate the proportion of 22 immune cells in GSE6613 and GSE98793, as shown in the graphs in Fig. 12A, B. The correlation of 22 types of immune cells revealed that CD8 T cells were negatively associated with Neutrophils (r = − 0.69) and that naive B cells were negatively related to plasma cells (r = − 0.61), whereas CD8 T cells were negatively related to naive CD4 T cells (r = − 0.57) in GSE6613 dataset. In addition, CD8 T cells were negatively associated with Neutrophils (r = − 0.67), and naive B cells were negatively related to memory B cells (r = − 0.69), whereas CD8 T cells were negatively related to naive CD4 T cells (r = − 0.61) in GSE98793 dataset (Fig. 12C, D). The immune cell infiltration of two datasets was compared in boxplots (Fig. 12E, F). The results showed that in the PD groups, only follicular helper T cells were lower than in the control group in the PD dataset. However, naïve B cells, plasma cells, CD8 T cells, and the other five cells were significantly different between the MDD groups and the control groups.
Besides, it was found that all hub genes were significantly correlated with neutrophil cells with p < 0.05 in GSE6613 (AQP9, r = 0.66, p = 2.  Scale bar = 500 μm. n = 6. B Represents data on immobility, climbing, and swimming times in the forced swimming test Interestingly, all hub genes exhibited a significant correlation with monocytes with p < 0.05 in GSE6613 (AQP9, r = 0.30, p = 5.0e-4; RPH3A, r = 0.17, p = 0.06; SPI1, r = 0.36, p = 3.0e-5, Additional Fig. 5G-I). These results suggested that the neutrophil cells and monocytes may play important roles in PD with MDD.

Discussion
As a progressive neurodegenerative disease, Parkinson's disease (PD) is characterized by bradykinesia, tremors at rest, rigidity, and postural instability. Depression, insomnia, constipation, and hyposmia, however, are not considered motor symptoms and often precede them [36]. The one of most common psychiatric symptoms of PD is depression (major depressive disorder, MDD), which has not been paid high attention to [37]. A study of meta-analyses found that MDD occurs in around 23% of patients with PD, which is higher than that seen in other chronic and disabling conditions [8,38]. The incidence of MDD will increase exponentially as PD progress [39]. MDD is rarely recognized as a manifestation, which delays diagnosis and treatment, and ultimately results in poor living conditions for patients. However, the associations and biomarkers of PD with MDD have not yet been clarified.
In the present study, we used bioinformatics analysis to construct a nomogram to assess MDD in PD patients and evaluated its diagnostic value. Public datasets of peripheral whole blood from PD and MDD patients are an efficient and practical method for clinical use. In addition to identifying three pivotal candidate genes (AQP9, RPH3A, and SPI1), we developed a nomogram for diagnosing MDD in PD patients.
Aquaporins 9 (AQP9), a novel biomarker found in our study, was used in the diagnosis of MDD in patients with PD. It belongs to the family of membrane proteins that mediate the transport of water [40]. AQP9 is the aquaporin most closely related to GlpF, a bacterial aquaglyceroporin that fluxes water as well as glycerol [41]. AQP9 is the first mitochondrial protein with two novel isoforms in brain tissue that is selectively expressed in dopaminergic neurons [42]. Studies in vitro showed that AQP4 and AQP9 transcription and expression changed dynamically after differentiation toward DAergic neurons, resulting in the vulnerability of the human SH-SY5Y cell line [43]. Additionally, in vivo, study demonstrated that the deletion of the aquaglyceroporin AQP9 is protective in a mouse model of Parkinson's disease [44]. Besides, AQP9 is also involved in the pathogen of epilepsy [45]. Above all, AQP9 plays an important role in the central nervous system including PD.
Rabphilin-3A (RPH3A) is a synaptic protein initially known as a synaptic vesicle-associated protein involved in the regulation of exo-and endocytosis processes at presynaptic sites [46]. RPH3A retains NMDA receptors at synaptic sites through interaction with GluN2A/PSD-95 complex [47]. The previous study has suggested that interfering with the formation of the RPH3A/GluN2A/ PSD-95 complex could lead to a more specific and direct approach to tackling the aberrant NMDAR localization and function in levodopa (L-DOPA)-induced dyskinesias [48]. Interestingly, a recent study has demonstrated that RPH3A overexpression improves motor behavior in the PD mice model, and is a novel molecular target to counteract α-syn-induced synaptic failure [49]. Besides, a study has suggested that the increase in Rph3A in the brain penumbra may be an endogenous protective mechanism against ischemia-reperfusion injury, which is mainly dominated by astrocytes [50]. RPH3A loss correlated with dementia severity, cholinergic deafferentation, and increased β-amyloid (Aβ) concentrations [51]. Additionally, RPH3A has also been reported that it associated with cognitive resilience via a proteome-wide association study of the human dorsolateral prefrontal cortex [52]. Briefly, RPH3A was involved in multiple neuropsychiatric diseases, which may as a target for PD with MDD.
SPI-1 proto-oncogene (SPI1, also known as PU.1) was found by Moreau-Gcherin et al. known as one of the erythrocyte transformation-specific (Ets) family of transcription factors (TF). A meta-analysis of the PD public dataset showed that SPI1 is the TF of the hub gene [53]. Meanwhile, SPI1 is also an important TF of the hub gene in PD with bipolar disorder [54]. Besides, a study has demonstrated that SPI1 is a key gene associated with postpartum depression (PPD) [55]. In addition, abnormalities in SPI1 in the brainspleen axis might, in part, play a role in the pathophysiology of MDD [56]. Therefore, SPI may as a key gene in PD with MDD. Notably, the database projections on transcription exhibit conformity with the qPCR outcomes, albeit exhibiting minor discrepancies at the protein level. More experimental verification should be given in the future.
One limitation of the present study was the lack of the common hub genes of PD and MDD that have been validated in patients with either disease, but not both. As PD and MDD patients are not included in the current datasets, a more comprehensive validation can be attempted in the future.

Conclusion
AQP9, RPH3A, and SPI1 have been identified as common hub genes in PD and MDD in this study, which was increased compared to the sham groups on transcription Fig. 11 PPI network and TFgene interaction network. A The gene-gene interaction network for hub genes was analyzed using the GeneMANIA database. The 20 most frequently changed neighboring genes are shown. B Red nodes represent high-confidence candidate genes, and green nodes represent TF genes in a network for TF-gene interaction level. However, AQP9 and SPI1 also increased while RPH3A was reduced via western blotting. Neutrophil and monocyte infiltration might play central roles in the development of PD and MDD, suggesting they may be potential targets for diagnosis and treatment.
Authors Contribution HQW analyzed the data and drafted manuscripts. SSD, CMW, and WMG participated in the revision of the manuscript and figures. BHC and FLY designed a significant research topic and revised the manuscript. All authors have read and agreed to the published version of the manuscript.

Funding
The present work was supported by the National Nature Science Foundation of China (No. 82071325) and the Jiangsu Province Capability Improvement Project through Science, Technology, and Education (No. ZDXK202215).

Data Availability
The data that support the findings of this study are available from public databases.