Identification of Gene Co-Expression Modules and Core Genes Related to Immune Disorders in Major Depression Disorder

Introduction

Depression is a psychological and mental illness that is very harmful to humans, often coinciding with an immune function disorder.¹ A number of patients receiving cytokine immunotherapy have been clinically observed to exhibit a series of depressive symptoms, and antidepressants can block these symptoms.^2–4 Animal studies have also revealed a series of changes in neuropsychological behavioral changes following the injection of inflammatory cytokines or endotoxins: anorexia, somnolence, lack of pleasure, weight loss, decreased activity, and even lack of cognitive function, commonly referred to as “pathological behavior”.^5,6 These pathological behavior and depressive symptoms are very similar to some extent, and their mechanisms are also similar with the mechanisms of stress-induced depression. Antidepressant effects are similar with the mechanism of pathological behavior, and antidepressants have a significant effect on pathological behavior. On the other hand, patients with chronic immune diseases often also present with depression, while anti-cytokines are associated with the improvement of depression accompanied by chronic immune diseases.^7–9 In 1999, the “inflammatory response model of depression” was proposed,¹⁰ which was considered to be associated with the activation of the inflammatory response system; it is a kind of psychoneuroimmunology disorder. The activation of the peripheral immune can lead to several behaviors related to depression as well as neuroendocrine and neurobiochemical changes by releasing inflammatory cytokines. The internal molecular mechanism of immune disorders and depression is not yet clear. The purpose of this research was to analyze the molecular mechanism and the key target of immune disorders in patients with depression by means of the combination of an immune scoring algorithm and WGCNA.

Methods and Materials

Data Collection

A data chip named GSE98793 was downloaded from the GEO database.¹¹ The study investigated the changes in gene expression in whole blood samples from donors diagnosed with severe depression (MDD) and from healthy controls. We could obtain microarray data of whole blood for the MDD patients (N = 128, 64 patients with extensive anxiety disorder; the diagnosis was made through a MINI questionnaire; 64 patients without anxiety disorder) and the healthy control group (N =64). The RMA algorithm was used to normalize the chip data to obtain the standardized gene expression values.

Evaluating the Proportion of Peripheral Blood Immune Cells

The proportion of different peripheral blood immune cells of patients with depression and normal volunteers was calculated by the Xcell database,¹² and then the proportions of different peripheral blood immune cells of normal volunteers and patients with severe depression were compared and analyzed, which was used as clinical data in co-expression network analysis.

Co-Expression Network Construction

Initially, all the samples and genes in the expression profile excelled at the test of high-quality genes or good samples. Then, a scale-free network of all genes was constructed by utilizing the WGCNA R software package.¹³ First, the genes with expression indicated as “NA” were removed according to the expression profile, the variance in each gene in each sample was calculated, the genes with a standard deviation greater than 0.2 were screened, all the samples were further analyzed by cluster analysis, and the outlier samples were eliminated according to the clustering distance. The scale-free co-expression network was realized when the soft threshold corresponded to 6 and the correlation coefficient increased to 0.8. The expression matrix was converted into an adjacent matrix, which was later transformed into a topological matrix. Premised on topological overlap matrix (TOM), we utilized the average-linkage hierarchical clustering technique to cluster genes, in accordance with the criterion for a mixed dynamic shear tree and established the minimum number of bases per gene network module at 30. After determining the gene module utilizing the dynamic shear method, we computed the eigenvector value for each of the modules in turn. Subsequently, we performed cluster analysis on the modules, combined close modules into new modules, and fixed height=0.25.

Enrichment Analysis of Key Co-Expression Modules

Gene ontology (GO) enrichment analysis was carried out utilizing the Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/home.jsp), which is thought to be the commonly utilized functional annotation tool.^14,15 The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway was employed as a reference, and the enrichment analysis was performed by DAVID. The genes in key co-expression modules with a high correlation with neutrophils were uploaded to study potential roles, and P < 0.05 was considered to represent significant differences.

Identification of Key Genes in the Modules

In the process of WGCNA of key genes in the modules, the interaction network of genes in the key module was derived, and the protein-protein interaction network was drawn by Cytoscape. The degree value of the key nodes in the network was calculated by Cytoscape,¹⁶ and the key nodes in the network were screened according to the degree value. In addition, the correlation between all genes in the module and the module could be analyzed, and the key genes in the module could be assessed in accordance with the correlation score between the genes and the module.

Results

Patients with Major Depression Disorder Show Visible Immune Cell Disorders

We utilized the Xcell database to analyze the scores of 11 common peripheral blood immune cells of patients with depression and normal volunteers. As shown in Figure 1, the proportion of CD4⁺ T cells, CD8⁺ T cells, B cells, and NK cells in the peripheral blood of MDD patients decreased significantly, while the proportion of macrophages in the peripheral blood of MDD patients increased. Previously, various studies have reported that depression patients tend to be immunosuppressed.^17–19 However, the molecular mechanism under immunosuppressed conditions is still unknown. It was reported that changes in the leukocyte subgroup in the peripheral blood of patients with depression, usually accompanied by disruption of neutrophil phagocytosis, inhibition of cytotoxicity of NK cells, and decrease in proliferation of lymphocytes. This is consistent with our analysis results. Furthermore, the increase in macrophage score may also be related to the increase in inflammatory factor level.^20,21 Overall, the characteristics of peripheral blood immune cells in depressed patients tend to be immunosuppressed, and the molecular mechanism needs to be further explored.

Figure 1 Comparison of immune cell fraction between normal volunteers and patients with major depression disorder.

Construction of a Peripheral Blood Gene Co-Expression Network for MDD Patients

To systematically explain the relationship between MDD patients and immune cells and the underlying mechanism. Firstly, we obtained expression profile data containing 128 samples and 21,138 genes, which was analyzed by cluster analysis. The clustering results with sample characteristics are illustrated in Figure 2A. The red representation in the figure is marked as a nonzero sample in the table Finally, by eliminating outliers, we obtained a novel data expression profile, including a total of 125 samples and 20,011 genes. We employed the R software package WGCNA to build a weighted co-expression network. The soft threshold of β = 6 was chosen as illustrated in Figure 2B in order to make sure that the network was scale-free. An aggregate of 13 modules was revealed by clustering. As shown in Figure 2C, these modules display and contain different numbers of genes in different colors, whereby, the gray module comprises a set of genes that could not be merged into other modules. Based on the eigenvectors of each module, we computed the correlation between these modules. As shown in Figure 2D, each module is separated from the other in the network, while there is a correlation between some modules. These results suggest that some modules have significant differences between each other, and some modules have similar expression patterns. The above data suggests that we have successfully constructed a gene co-expression network of blood samples from MDD patients, which can be used for subsequent in-depth analysis.

Figure 2 Construction of a WGCNA co-expression network of MDD patients. (A) Sample clustering molecular map; red representation of samples is marked as nonzero in the table. (B) Calculation and selection of the optimal soft threshold power. The influence of different powers on scale Independence (left), and the effect of different powers on mean connectivity (right). (C) Tree diagram of clustering of 20,011 genes and 13 combined modules in 125 MDD patients. (D) Correlation heat map of 13 modules.

Modules Related to Immunosuppressive Status in MDD Patients

In order to further understand which genes can affect the peripheral blood immunosuppression of depressed patients, we computed the relationship between clustered modules and each phenotype index (Immune cell scores). As depicted in Figure 3, we can see the correlation and statistical significance between the phenotype and each module. Among them, CD4⁺ T cells, CD8⁺ T cells, B cells, and NK cells were positively correlated with the yellow module but negatively correlated with the purple module. Immune cells scores in the peripheral blood of MDD patients were considerably reduced than those in the healthy volunteers. In addition, we also analyzed the correlation of the module genes and phenotypes. As shown in Figure 4, we can see the correlation between the CD4⁺ T cells, CD8⁺ T cells, B cells, and NK cells with each module, in which the yellow and purple modules are the two modules with high correlation. Figure 5 shows the correlation scatter diagram of several immune cells (CD4⁺ T cells, CD8⁺ T cells, B cells, and NK cells) with genes in interested modules, which show a positively correlated trend. All in all, we have explored the connection between key immune cells and key genes from different angles, in order to further lock key genes.

Figure 3 Correlation heatmap of Co-expression module and immunocyte score of MDD patients. Heatmap of the correlation between each module and the blood immune cell scores of MDD patients, red represents positive correlation, blue represents negative correlation. Each unit contains the corresponding correlation coefficient and P-value.

Figure 4 The correlation histogram of genes in the modules and immunocyte scores. There were significant levels of 13 co-expression modules related to B cells, CD4+ T cells, CD8+ T cells, and NK cells.

Figure 5 Scatter diagram of the correlation between genes and immunocyte scores. Scatter diagram of the correlation between genes and four immunocyte scores (B cells, CD4+ T cells, CD8+ T cells, NK cells) in the yellow and purple module.

Gene Function Enrichment Analysis in Key Modules Related to Immune Disorder

As mentioned above, the specific molecular mechanism of immunosuppression in depressed patients is not yet clear, so we conducted an in-depth analysis of the yellow and purple models to explore their potential molecular mechanism. According to DAVID enrichment analysis of genes in the yellow and purple models, as illustrated in Figure 6A, the yellow module is predominantly concentrated in thymic development, T cell co-stimulation and differentiation, B cell activation, and the peripheral nervous system myelin maintenance. The KEGG signaling pathway analysis for the yellow module is mainly enriched in primary immunodeficiency signaling pathways, as seen in Figure 6B. While, genes in the purple module are mainly responsible for protein phosphorylation, leukocyte migration, and promoting apoptosis (Figure 6C). They also relate to the Rap1 signaling pathway (Figure 6D) and the HIF-1 signaling pathway. Some studies have shown that the HIF-1 signaling pathway takes part in immunosuppression.²² The above analysis shows that genes in the yellow module are more connected to the positive regulation process of immunity, while those in the purple module are mainly related to the process and pathway of immune negative regulation. The following work can focus on the core genes in these two modules, and the study of these genes may provide new ideas for the treatment of MDD patients.

Figure 6 GO and KEGG pathway enrichment for genes in the immune related modules. (A) The top 15 enriched biological processes of genes in the yellow module. (B) The top 5 enriched signaling pathways of genes in the yellow module. (C) The top 15 enriched biological processes of genes in the purple module. (D) The top 5 enriched signaling pathways of genes in the purple module.

Key Genes Related to the Immunosuppressive State in MDD Patients

To discover the core genes in the yellow and purple modules, we performed protein network analysis on these two modules respectively. The PPI network of the yellow module contains 583 genes and 49,683 interaction relations (Figure 7A). The purple module contains 79 genes and 1827 interaction relations (Figure 7B). As shown in Figure 7A and B, nodes’ colors are connected to the degree value of each node. The red color denotes a higher degree while the green color denotes a lower degree. Finally, we defined the 10 genes with the uppermost degree value in the network as key genes as previously described.²³ In addition, module genes with MS score greater than 0.85 were also included as key nodes. As shown in Figure 8, the intersection of the two criteria was together used to find the key genes in each module. The key genes in the yellow module are IARS, SKAP1, LCK, PEBP1, FAM162A, ILF3, LDHB, RPSA, CD247, and EEF1B2. LDHB encodes the B subunit of the lactate dehydrogenase enzyme and is vital in the anaerobic conversion of lactate to pyruvate, NADH to NAD^+, and conversion between NAD⁺ and NADH. An increase in LDHB protein expression has been identified in peripheral blood mononuclear cells of initial onset antipsychotic naive schizophrenics.²⁴ GWAS analysis based on MDD patient samples showed that significant GWS signals were found around the LDHB gene in European American samples, suggesting that there is a close association between LDHB and MDD.²⁵ Furthermore, SKAP1 is an adaptor with a peculiar N terminus, a C terminal SH3 domain, and a PH domain.²⁶ SKAP1 can act as a scaffold for PLK1 in regulating the cell cycle of T-cells.²⁷ Existing studies indicate that SKAP1 has a positive regulatory effect on immunity.^28,29 In addition, LCK, RPSA, and CD274 have been reported to be related to immune activation.

Figure 7 Hub genes in immune related key modules. (A) Gene PPI analysis, visualization, and 10 core genes in the yellow module. (B) Gene PPI analysis, visualization, and 10 core genes in the purple module.

Figure 8 Core genes in immune related key modules. The key genes in the PPI network and the genes whose correlations with modules were higher than 0.85 were analyzed by a Venn diagram and screening of the core genes that were consistent in both. The top is the Venn diagram of the yellow module, and the bottom is the Venn diagram of the purple module.

The most critical genes in the purple module are LOC101928143, LYN, RTN3, DOCK5, RALB, REPS2, ACOX1, CFLAR. In a study involving patients with MDD, SNP analysis showed a rare heterozygous mutation in exon 31 of DOCK5 (c.3170A>G, p.E1057G). The convergent functional genomics analysis illustrated that DOCK5 was recorded among the biomarkers for suicidality and mood state.³⁰ The RAS-like GTPase RALB employs SEC5 to trigger innate immunity signaling and build a bridge between inflammation and immunity.^31,32 Furthermore, LOC101928143, in the Purple module is an RNA gene, which belongs to the lncRNA class.³³ Currently, there is a lack of reports on the function of this gene, which deserves further study. In general, we have screened some genes that are closely related to MDD, some of which play a role in immune activation, and some play a role in immunosuppression, which provides the possibility for subsequent in-depth research.

Conclusion

This study confirmed that the proportion of CD4⁺ T cells, CD8⁺ T cells, B cells, and NK cells in the peripheral blood of patients with depression were noticeably lower than that in the peripheral blood of normal volunteers by the immune scoring algorithm. As depression is the comorbidity of many diseases including thyroid disease,³⁴ cancer,³⁵ diabetes,³⁶ cardiovascular³⁷ and nervous system disease. Immune cells such CD4⁺ T cells, CD8⁺ T cells, B cells, and NK cells also play an important role in other disease progression. The pathogenesis of polycystic ovary syndrome is related to immune cells including T, CD4⁺ T, and NK cells.³⁸ Moreover, patients with pancreatic ductal adenocarcinoma M0 macrophages, gamma delta T cells and naive CD4 T cells were independent predictive factors of a poor outcome for PDAC patients.³⁹ Similar to these diseases, our research also confirmed that these types of immune cells play an important role in the progression of MDD. In addition, we found gene expression modules closely related to the decrease of these cells by WGCNA and screened and analyzed the key genes in these two modules. Among them, SKAP1 and RALB may be essential regulatory genes affecting the immune status of patients with depression, which is worthy of further experimental study.