Analysis of Differentially Expressed Genes That Aggravate Metabolic Diseases in Depression

Depression is considered the second leading cause of the global health burden after cancer. It is recognized as the most common physiological disorder. It affects about 350 million people worldwide to a serious degree. The onset of depression, inadequate food intake, abnormal glycemic control and cognitive impairment have strong associations with various metabolic disorders which are mediated through alterations in diet and physical activities. The regulatory key factors among metabolic diseases and depression are poorly understood. To understand the molecular mechanisms of the dysregulation of genes affected in depressive disorder, we employed an analytical, quantitative framework for depression and related metabolic diseases. In this study, we examined datasets containing patients with depression, obesity, diabetes and NASH. After normalizing batch effects to minimize the heterogeneity of all the datasets, we found differentially expressed genes (DEGs) common to all the datasets. We identified significantly associated enrichment pathways, ontology pathways, protein–protein cluster networks and gene–disease associations among the co-expressed genes co-expressed in depression and the metabolic disorders. Our study suggested potentially active signaling pathways and co-expressed gene sets which may play key roles in crosstalk between metabolic diseases and depression.


Introduction
Major depressive disorder (MDD) is a psychiatric disease. It is considered one of the major global health burdens [1]. An estimation by the WHO suggested that 350 million people of all ages suffer from MDD that involves serious suicidal ideation [2,3]. External stressors can act as stimuli that cause MDD, depending on their repetitiveness and durations [4]. The pathophysiological responses in MDD include adaptive hormonal changes, such as increases in corticosteroids and adrenocorticotropin, which altogether play a significant role in the hyperactivity of the hypothalamic-pituitary-adrenal (HPA) axis at the onset of depression [5]. Studies over the years have found that there is a significant association between obesity and depression [6]. The excess intake of macronutrients above physiological requirements raises the sugar and the fat depots by increasing body weight, which promotes obesity. Obesity has been associated with many metabolic dysfunctions which in turn may act as stimulators of MDD [7]. Metabolic disorders such as hyperlipidemia and hyperglycemia are prevalent in a wide range of the psychiatric patients with MDD [8]. However, a strong link between depression and obesity has not been established [6]. Reports suggest that prediabetic patients and undiagnosed diabetic patients have more prevalence of depression compared to nondiabetic individuals [9]. Patients with type 1 and type 2 diabetes tend to acquire much more clinical symptoms of depression than nondiabetic individuals [10]. Overactivation of innate immunity by the cytokine-mediated inflammatory response leads to dysregulation of the hypothalamic-pituitary-adrenal (HPA) axis [11]. Hyperactivation of the HPA axis can disrupt metabolic homeostasis and promote obesity. HPA axis hyperactivation can include corticotrphin releasing hormone (CRH) suppression and leptin resistance. Activation of the HPA axis and the sympathetic nervous system (SNS) increases the production of adrenalin, noradrenalin and cortisol, promoting insulin resistance, visceral obesity and diabetes. Simultaneously, proinflammatory cytokines might directly affect the brain, causing depressive symptoms [12]. Reports from meta-analyses also support that depression may increase the risk of developing type 2 diabetes [13]. Hyperphagia is a commonly observed symptom in patients with depressive disorders, and this may cause issues with hepatic biochemistry, including hepatic injury, elevated hepatic enzymes and loss of hepatic blood flow [14]. These histological abnormalities during depression may cause liver disease: NAFLD ranges from simple steatosis to NASH and fibrosis [15]. Reports showed that the prevalence of MDD is higher in the NASH population than that in control subjects [16]. Importantly, they found that the diagnosis of MDD tended to be associated with steatosis grades [17]. It is widely accepted that psychiatric disorders, including MDD, are associated with histological severity of liver diseases [17]. As it is known that depression is associated with inflammation [18], which in turn affects insulin resistance, which could be relevant to the causality of NASH [14]. The associations between depressive disorder and the above-mentioned metabolic diseases have been established, but the genetic co-expression and interconnections among them have not been well elucidated yet. We employed a systematic approach to study independent microarray datasets related to depression, obesity and related metabolic disorders-diabetes and NASH. To understand the effects of depression on these metabolic diseases, we studied differentially dysregulated genes, enriched pathways, ontology analysis and protein-protein interactions; and lastly, for validation of our study, we compared the results with another OMIM disease database and datasets related to depression.

Overview of the Study
To identify links between depression and its effects on metabolic disorders, we employed microarray datasets from NCBI. Differentially expressed genes (DEGs) identified in these datasets were further analyzed to find out the commonly expressed DEGs among depression and related metabolic diseases. These genes were further used to construct a gene-disease association network, analyze KEGG pathways, analyze ontological pathways (GO) and perform protein-protein-interaction (PPI) network analysis. To validate the potential DEGs identified from enrichment analysis, the OMIM disease database and dbGap were used.

Retrieval of Microarray Datasets
We analyzed gene expression data with 4 independent experimental microarray datasets-accession numbers GSE58430 [19], GSE128021 [20], GSE43950 [21] and GSE43600 [22]-from the National Center for Biotechnology Information (NCBI) GEO2R database, as shown in Table 1. Dataset GSE58430 was produced by using human transcriptomic profiling of peripheral blood CD4+ T-lymphocyte cells. There were four groups: nondepressive asthma, depressive asthma, depression and healthy controls. We compared depressed patients with the healthy controls. The obesity dataset GSE128021 was generated from gene expression differences in omental mesothelial cells from the lean and obese human donors with a cross-sectional case-control study with two different cohorts of lean and obese patients with varying different degrees of obesity. The diabetes dataset GSE43950 was produced by gene expression profiling in endothelial precursor cells of patients protected from microvascular complications by modulating the TGF-β/PAI-1 axis in CD34+ cells from diabetic patients and controls. The sample of the dataset was constructed with peripheral blood from healthy controls, and from patients with longstanding, poorly controlled diabetes with severe microvascular complications blood and CD34+ cells were obtained. RNA extraction was followed by AffyNugen amplification, and the cDNA was probed to the Human RSTA Affymetrix 2.0 chip. The NASH dataset Life 2021, 11, 1203 3 of 17 GSE43600 was constructed with liver samples from patients with varying severities of steatosis. The summarized description of the datasets is shown in Table 1.

Differential Gene Expression Analysis and Validation with Random Unrelated Diseases
Each dataset was normalized by quantile normalization, and the R package limma was used to identify the DEGs between depression and related metabolic disease datasets. The batch effects from each dataset were removed by ComBat method [23]. To determine the upregulated and downregulated genes, the threshold values were set to log FC > 1 p value < 0.05. The Venn diagram was also constructed using the Venn Diagram package in R. All significant DEGs are presented in volcano plots generated using R software.
To validate the datasets adopted, we have compared selected diseases in the present study with an unrelated, random infectious disease, influenza, shown in Table 2. Interestingly, we found that below 2% of the upregulated and 1% of the downregulated genes were observed in the datasets for unrelated diseases. Altogether, these findings suggest that the observed common DEGs were not found by chance, but the relatedness of the diseases caused higher chances of their occurrence; hence the rationale of this study is to find the DEGs expressed commonly between depression and metabolic disorders.

Gene Set Enrichment Analysis
Pathway-based analysis reveals the molecular crosstalk in the progression of complex diseases networks. We analyzed pathways of the commonly altered DEGs in depression and related metabolic diseases datasets using Enrichr [24], a comprehensive web-based gene-set enrichment tool, to construct KEGG pathways. Enrich R implements the Fisher exact test, in which a binomial distribution and independence are assumed for the probability of any gene belonging to any set from the random input gene list, in order to create the mean rank and standard deviation from the expected rank. It is followed by the correction to the next Fisher exact test, which calculates a Z score for the standard deviation. This is considered as a new, corrected score. Alternatively, a p-value obtained from the Fisher exact test can be combined with the Z score of the deviation. Hence, we combined the p-value obtained from the Fisher exact test and combined with the z-score of the deviation from the expected rank by multiplying these two numbers as follows: c = log(p)·zc = log(p)·z where c is the combined score, p is the p-value computed using the Fisher exact test, and z is the z-score computed by assessing the deviation from the expected rank.
The clustering level z-scores and p-values are highlighted in red if the clustering is significant (p-value < 0.1). This clustering indicator provides an additional assessment of related genes and measures their relevance in the specific gene-set libraries from the input list of genes. The observation of one or two clusters on the grid suggests that a gene-set library is relevant to the input list. It also indicates that the terms in the clusters are relevant to the input list. We considered significant signaling pathways after applying several statistical analyses in R packages, such as fgsea and MSigDB to access the KEGG gene sets, and clusterprofiler where we used hypergeometric tests to find out the statistically enriched KEGG pathways. To prevent a high false-discovery rate (FDR < 0.05) in multiple testing, q values were also estimated. We performed gene-set enrichment analysis for KEGG pathways significantly associated with upregulated and downregulated DEGs which were identified from the employed datasets (p value < 0.05)

Ontology Pathway
To analyze the ontology pathway significantly, potential biological process (BP), molecular function (MF) and cellular components (CC) involved in overlapping DEGs among depression and other metabolic diseases, we used the online database Enrichr to conduct the ontology-pathway enrichment analysis.

Protein-Protein Interaction (PPI) Network and Hub-Gene Identification
The protein-protein interactions of overlapping DEGs with a combined score > 0.4 were identified with systematic approaches in the STRING database [25], and we mapped the protein functional associations and protein-protein interactions (PPI). We use Cytoscape software [26] (http://www.cytoscape.org/, (accessed on 20 June 2021), version 3.7.1; Institute for Systems Biology, Seattle, WA, USA software plugin) to visualize and construct the transcriptional regulatory network of common DEGs among the metabolic disease and depression database. After overlaying DEGs on networks in CytoHubba plugin in Cytoscape, hub genes with >10 degree were identified among all the modules from different datasets.
Molecular complex detection was used to find the most significant functional interactions between proteins' modular clusters from dense PPI network regions [27]. Module identification criteria were included with a degree cut-off of 2, node score cut-off of 0.2, k-core of 2 and the maximum depth of 100. Significant modules were identified with MCODE score > 3 and nodes > 3. The cluster networks were visualized with Cytoscape 3.7.1 software plug-ins.

Prediction of the Master Transcription Factors (TFs)
To predict master TFs that significantly regulate the DEGs, we have utilized the iRegulon plugin of Cytoscape software (version 3.8.0) to detect regulons from all DEGs [28]. The iRegulon method is determined by a ranking-and-recovery system. In this method, all genes of the human genome are scored by a motif discovery step which connects the cluster binding sites within cis-regulatory modules (CRMs) and the potential distal location of CRMs of the transcription start site (TSS ± 10 kb) in both upstream and downstream. The normalized enrichment score (NES) is computed at the recovery step for TFs of each set of genes. The prediction of the TFs is based on NES and their putative target genes directly from the input lists. The association of TFs to motifs using both explicit annotations and predictions of TF orthologs and motif similarity is optimized by this method. A transcription factor normalized enrichment score was computed for each group, where a normalized enrichment score > 4.0 is considered significant, and the maximum falsediscovery rate (FDR) for motif similarity was set as 0.001.

OMIM Disease and OMIM Expanded Database
In order to validate the results obtained from various analyses, we selected the OMIM database [29] and the dbGap benchmark database [30] to verify the significantly expressed genes from co-expression analysis in depression and other metabolic disease datasets.

Validation of the Common DEGs with other Depression-Related Datasets
To compare the DEGs commonly expressed on each disease-related dataset, we compared the expression in other datasets related to depression and depression-induced carcinoma. GSE 14922 [31] is produced by transcriptional profiling of human adrenocortical tumors in which we compared primary adrenocortical cancer tissues with the normal adrenal cortex. GSE109857 [32] is composed of WHO-classified grade-III glioma with the healthy controls. In GSE114852 [33], prenatal exposure to maternal stress and depression has been considered as the risk factor for depression, and we compared healthy mothers with those with depression. In GSE32280 [33], we compared the whole-genome microarray expression profile of leucocytes with SSD (somatic symptom disorder) and MDD patients.

Statistical Analysis
Statistical analysis was performed using the Origin 8 software (version 8.6, OriginLab Corporation, Massachusetts, MA, USA) and Graphpad Prism 8 (GraphPad, San Diego, CA, USA). The unpaired Student's t-test (two-tailed) and one-way ANOVA followed by Fisher's least significant difference (LSD) post hoc test were performed to analyze the data, where appropriate. A value of p < 0.05 was considered as statistically significant. N.S.: not significant (p > 0.05).
To avoid complications arising from the different experimental systems, compound traits or ethnicity used in the original datasets, we normalized the gene expression data by using Z-score transformation (Zij) for each type of tissue gene expression profile Zij = (gij− mean g)/(SD(gi)) where SD represents the standard deviation, gij denotes the value of the gene expression i in sample j. Transformed gene expression values using the Z score of different diseases at different platforms can be compared according to a previous study [34].

Identification of DEGs
We studied selected, independent, microarray datasets, and after normalization of the dataset GSE58430 (depression) showed 1092 DEGs (345 up,747 down), the dataset for obesity, GSE128021, contained 2422 DEGs (1190 up and 1230 down), the dataset for diabetes, GSE43950, showed 1179 DEGs (848 up, 331 down), the dataset for NASH, GSE46300, showed 1577 DEGs (1227 up, 350 down), which are listed in Supplementary Data S1-S4, respectively. Significantly expressed upregulated and downregulated DEGs for individual datasets were plotted in volcano plots as shown in Figure 1. In Figure 2, the Venn diagram made in R represents the number of co-expressed DEGs among the individual datasets. We found that depression datasets share a number of differentially expressed genes with the datasets of obesity (113 genes), diabetes (83 genes) and NASH (71 genes). We tried to emphasize the co-expressed genes among the associated disease between depression and metabolic disease groups. Comprehensive bioinformatics was performed to determine the independent DEGs across all the datasets. Enrichment pathway analysis was also performed to identify the pathways between the co-expressed DEGs. Notably we found five genes to be dysregulated among all the selected datasets namely, PRDM2, CXCL1, PHLDA1, DIDO1, CDA. Figure 3A-C represents the heatmaps of commonly altered genes among depression and obesity, diabetes and NASH, respectively to identify the expression pattern according to their changes in log fold values. with depression and three other metabolic disease-related datasets and found the significantly enriched ontology pathways in the DAVID database [35]. Tables 3-5 showed the ten most significant ontology pathways linked with biological pathways, molecular functions and cellular components, respectively, in obesity, diabetes and NASH.

KEGG Pathway Analysis
Enriched pathways are a series of molecular mechanisms and their interconnections. To identify a disease network pathway that plays an important role we employed commonly altered DEGs in depression and three other metabolic diseases to identify the biologically active pathways. After applying several statistical analyses, the ten most significantly enriched pathways in KEGG are shown in Figure 4. KEGG pathways associated with depression and diabetes or NASH are mostly related to infectious diseases including legionellosis, hepatitis C tuberculosis, measles, and salmonellosis. Notably, some of the pathways affect TNF and IL-17 signaling pathways, phospholipase D signaling pathways, protein digestion and absorption, and glycophospholipid metabolism. Among all the significantly enriched pathways, the TNF signaling pathway is commonly found in obesity and diabetes, whereas the IL-17 signaling pathway and sulfur relay system are common in both obesity and NASH.

Gene Ontological Pathway Analysis
Ontology pathways includes all the altered biological and disease pathways and their relationship in the disease domain. We utilized the commonly altered genes associated with depression and three other metabolic disease-related datasets and found the significantly enriched ontology pathways in the DAVID database [35]. Tables 3-5 showed the ten most significant ontology pathways linked with biological pathways, molecular functions and cellular components, respectively, in obesity, diabetes and NASH.
The significant pathways associated with the common DEGs in depression and obesity are the type I interferon signaling pathway, cytokine-mediated signaling pathway, and insulin secretion involved in the cellular response to glucose stimulus. For common DEGs in depression and diabetes, we found interleukin-1 receptor binding, peptidoglycan binding, phosphatidylinositol-4,5-bisphosphate 3-kinase activity, phosphatidylinositol bisphosphate kinase activity, chemokine receptor activity, inflammatory response, and cytokine-mediated signaling pathway. For common DEGs in depression and NASH, we found the cardiolipin biosynthetic pathway, phosphatidylglycerol biosynthetic pathway, thyroid hormone pathway and cardiolipin metabolic pathway. Mostly, these pathways are associated with lipid biosynthetic pathways and inflammation pathways.

PPI Interactions
A protein-protein interaction helps to establish the connection between two specific proteins to build up functional biochemical networks for biological functions in the cells. We constructed protein-protein interaction networks with commonly altered genes associated with depression and the other three metabolic diseases with the STRING database [25], and visualized clusters consisting of hub genes discovered in each experimental dataset using the cytoscape software, as shown in Figure 5. In the CytoHubba plugin of cytoscape software altogether 29 hub genes were discovered, having >10 degrees and maximum nodes and neighborhood components. Among the 29 hub genes, 5 hub genes were differentially expressed across all the experimental datasets including depression, obesity, diabetes and NASH.

MCODE
Molecular complex detection was used to find the most significant functional interactions between proteins modular clusters from the dense PPI network region. PPI networks in obesity and depression datasets showed two clusters, as shown in Figure 6. Cluster 1 contains four nodes with PTGS2, PTX3, CDA, CXCL1 genes and cluster 2 contains three nodes with FUT4, MME and CD2 genes. Simultaneously, two clusters were obtained from the PPI interactions for depression and NASH datasets, and both contain only three nodes. One contains CUX1, STAT3, EOMES and the other one contains ACTA2, FLNA and WDR1. For diabetes and depression datasets there are three clusters with multiple nodes. Cluster 1 contains six nodes with CXCR1, IL1RN, CXCR2, IL1B, FPR1, CLEC4E, with 10 edges and four density nodes. Cluster 2 contains NOD2, CXCL1, QPCT, CDA, PTX3, TLR2 with ten edges and four density nodes. Cluster 3 has only three nodes, with XAF1, IFIT2, SAMD9L genes.

TFs Prediction by iRegulon
To identify gene regulatory networks we used iRegulon, a computational method used to discover various transcription factors and their target genes, which in turn can identify a set of co-expressed genes. In the common genes between depression and obesity, we found five TFs. Among them FOXN4, SPDEF, NKX2-1 have 12, 10 and 9 targets, respectively, and MECOM and HDAC2 have 7 targets. MECOM is identified as an oncogene, stimulating cell proliferation and development, whereas HDAC2 is an important gene that belongs to the histone acetylase family and plays a major part in cell-cycle progression and developmental events. Among the common genes of depression and diabetes, there are four TFs, namely, BCL3, MXI1, GMEB2, NFKB1. BCL3 and GMEB2 have 31 targets. BCL3 is a proto-oncogene, and it acts as a transcription coactivator and activates GMEB2 through NFKB dimers, which increases sensitivity to lower the concentrations of glucocorticoids. The common DEGs of depression and NASH have four TFs, i.e., ATF, SRF, SREBF1, SCRT2. SREBF1 and SRF have the highest targets with 18 and 14, respectively. SREBF1 promotes cholesterol biosynthesis and lipid homeostasis. On the other hand, the SRF transcription factor is identified as a novel upstream mediator of stress.

TFs Prediction by iRegulon
To identify gene regulatory networks we used iRegulon, a computational method used to discover various transcription factors and their target genes, which in turn can identify a set of co-expressed genes. In the common genes between depression and obesity, we found five TFs. Among them FOXN4, SPDEF, NKX2-1 have 12, 10 and 9 targets, respectively, and MECOM and HDAC2 have 7 targets. MECOM is identified as an oncogene, stimulating cell proliferation and development, whereas HDAC2 is an important gene that belongs to the histone acetylase family and plays a major part in cell-cycle progression and developmental events. Among the common genes of depression and diabetes, there are four TFs, namely, BCL3, MXI1, GMEB2, NFKB1. BCL3 and GMEB2 have 31 targets. BCL3 is a proto-oncogene, and it acts as a transcription coactivator and activates GMEB2 through NFKB dimers, which increases sensitivity to lower the concentrations of glucocorticoids. The common DEGs of Life 2021, 11, 1203 12 of 17 depression and NASH have four TFs, i.e., ATF, SRF, SREBF1, SCRT2. SREBF1 and SRF have the highest targets with 18 and 14, respectively. SREBF1 promotes cholesterol biosynthesis and lipid homeostasis. On the other hand, the SRF transcription factor is identified as a novel upstream mediator of stress.

Validation by OMIM Disease and dbGap
For validation of the DEGs obtained from PPI interactions and significantly enriched pathways in the depression dataset, we compared all the associated genes with the dbGaP, OMIM Disease and OMIM Expanded databases using only differentially expressed genes of depression so that we could potentiate the association between metabolic diseases and depression. After the analysis of all the related diseases, we found our selected three diseases, i.e., obesity, diabetes and NASH, among the list of all metabolic diseases which are associated with depression, as depicted in Figure 7.

Validation of the Commonly Expressed DEGs with Other Depression-Related Datasets
Our results show that PRDM2 is upregulated in two other datasets while it is downregulated in another. CXCL1 is downregulated in all other datasets including the dataset we analyzed. PHLDA1 is upregulated in all the datasets, which is inconsistent with our result. Both DIDO1 and CDA are upregulated in our dataset which is consistent with two other datasets, but downregulated in the other two, as referred to in Table 6.

Validation of the Commonly Expressed DEGs with Other Depression-Related Datasets
Our results show that PRDM2 is upregulated in two other datasets while it is downregulated in another. CXCL1 is downregulated in all other datasets including the dataset we analyzed. PHLDA1 is upregulated in all the datasets, which is inconsistent with our result. Both DIDO1 and CDA are upregulated in our dataset which is consistent with two other datasets, but downregulated in the other two, as referred to in Table 6.

Discussion
In our study, all the disease conditions selected are interconnected. We first hypothesized that the occurrence of one disease type may increase the chance of occurrence of the remaining three. The cellular origins of microarray datasets compared in our analysis are very diverse. Therefore, it cannot be ruled out that the expression values of a specific gene could be influenced by heterogeneous cell lines. The effect of these diseases is tissue-specific and dominant in those specific cell lines, hence, to mask the complications arising from the different experimental systems, compound traits or ethnicity employed in the original datasets, we normalized the gene expression data using Z-score transformation.
In this study, we investigated the effect of depression on its comorbidity in three metabolic diseases, i.e., obesity, diabetes and NASH. Hence, we compared the co-expression of the up-and down-regulated genes on MDD patients with that of obese, diabetic and NASH patients. To identify the commonly altered genes and their significantly enriched signaling pathways between depression and metabolic diseases, our study aimed to analyze the GEO microarray datasets on patients with depression, obesity, diabetes and NASH and their control datasets.
Our gene expression analysis suggested that depression and these metabolic diseases share a good number of commonly dysregulated genes. This indicates that depression has a strong influence on metabolic diseases [36]. The HPA axis is affected in depression, as a result it promotes hyperphagia which has a direct effect on metabolism. On the contrary, metabolic diseases cannot affect the HPA axis, so we cannot conclude that depression has a strong influence on metabolic disease or vice versa. To elucidate the molecular aspects of the regulatory network, we focused on pathway-based analysis, which is considered a new approach, to elucidate the molecular disease-network complexities related to the onset of depression. We found significantly enriched KEGG pathways for commonly dysregulated genes in depression related to each metabolic disease. These identified pathways strongly suggest that depression acts as a pivotal risk factor aggravating several metabolic disorders. Notably, gene expression ontologies and protein-protein interactions of commonly altered genes in depression and metabolic diseases revealed that depression plays a vital role in the prognosis of several metabolic diseases. We obtained two clusters among the common DEGs of both depression and obesity, and depression and diabetes, whereas three clusters were obtained in depression and NASH. PTX, CDA, and CXCL were present in the clusters among depression and obesity, and depression and NASH.
We verified our results with the gold benchmark database and found several genes play significant roles in both depression and other metabolic diseases. We collected disease names from the OMIM diseases, OMIM expanded, dbGap database using the commonly expressed dysregulated gene sets obtained from co-expression results in depression and metabolic diseases. Interestingly we identified the three metabolic diseases in the OMIM disease databases which coincide with the metabolic diseases we selected for the analyses. Moreover, the identified genes were significantly matched with the data found in these datasets.
Additionally, we investigated five genes, PRDM2, CXCL1, PHLDA1, DIDO1, CDA, which were commonly expressed in all datasets including depression, obesity, diabetes, and NASH. In our results, we found that PRDM2 is downregulated in NASH and depression datasets whereas it is upregulated in diabetes and obesity. On the other hand, except for obesity, all other datasets showed upregulation of CXCL1. PHLDA1 is upregulated in all the datasets except obesity. DIDO1 is only upregulated in the NASH dataset, and lastly, CDA is downregulated in obesity and NASH and upregulated in the other two datasets. Supplementary Data S5 shows the expression levels of these five commonly altered genes in each experimental dataset. To validate these commonly expressed DEGs among our experimental datasets, we compared their expression with other microarray GEO datasets related to depression in humans to ensure that the expression of these five genes were significantly related to the onset of depression. The expression level of these DEGs varied in all datasets related to depression due to the variation in the diverse cellular origins and the sample sizes. Our results are also in accordance with previously reported studies that suggest the same physiological role of these genes regarding depression and the metabolic diseases which we have obtained from our result. W. Fang et al. suggested that chromosomal deletion in the RIZ locus in PRDM2 may contribute to hepatocellular carcinoma which concomitantly can affect the progression of NASH [37,38]. The role of PRDM2 in high-fat induced obesity has been reported by Xie et al. via AKT transcription and AKT phosphorylation pathways [39]. Another chemokine ligand, CXCL1, plays a pivotal role in the pathogenesis of depressive disorder via inflammatory cytokines, as suggested by Hui et al. andŚlusarczyk et al. [40,41]. The CXCL1 chemokine gradient is also considered as a key factor to mediate obesity-dependent tumor-growth promotion according to the study by Zhang et al. [42]. In another study, CXCL1 was identified as one of the important markers of islets dysfunction and failure in T2DM [43]. CXCL1 is also elevated in NASH due to neutrophil infiltration [44]. Another important commonly dysregulated gene in our study, PHLDA1, is responsible for postpartum depression [45]. Loss of PHLDA1 also causes obesity, insulin resistance and NASH, while regulating lipogenesis, as suggested by Basseri et al. [46]. On the other hand, DIDO1 is associated with depression by affecting thyroid hormone levels [47]. Lastly, another co-expressed gene we found in our study was CDA, whose gene deficiency can lead to replicative stress as suggested by Farnces et al. [48].
It should be noted as a limitation of our study was that our report analyzed publicly available independent datasets which contain different types of cells and sample sizes. Consequently, it should also be considered that different cell lines have different expression values for a specific gene. Hence, this area should be further explored with in vivo and in vitro analysis for depression-related dysfunction, particularly with the brain region, which may strengthen the conclusion.

Conclusions
In this study, we considered gene expression (GEO) microarray data from depression and obesity, diabetes, NASH and control datasets to analyze and investigate the genetic links between depression and its effects on metabolic disorders. We analyzed gene expression, constructed gene-disease association networks, identified signaling and ontological pathways, analyzed protein-protein interaction networks, validated our results from the OMIM disease database, and finally compared the differentially co-expressed genes' expression with other GEO datasets. The outcome of our study confirmed that depression may exert a strong influence on metabolic disorders. Moreover, we identified five potential DEGs that are co-expressed in all the experimental datasets in depression and metabolic diseases. We believe our study will be useful for making more accurate disease predictions and identifying potentially better therapeutic approaches.

Data Availability Statement:
The data that support this study are available within the reference and its supplementary data files or available from the authors upon request.