Identification of Hub Genes in Comorbidity of Psoriasis and Vitiligo Using Bioinformatics Analysis

Introduction

Psoriasis is a chronic, immune-mediated inflammatory skin disorder.¹ The most widely recognized mechanism is the tumor necrosis factor-α/interleukin-23/interleukin-17 axis.² It is often characterized by widespread scaly erythematous plaques, which cause severe physical and psychological burden severely.^3,4 Patients with psoriasis are at an increased risk of developing other autoimmune diseases such as autoimmune bullous diseases, vitiligo, alopecia, and thyroiditis,⁵ and also chronic health problems such as cardiovascular diseases, diabetes, depression, and metabolic syndrome.^6–8

Vitiligo is an acquired depigmentation of skin mucosal disease.⁹ The mechanism underlying vitiligo is mainly related to melanocyte destruction or dysfunction. The frequency of concomitant autoimmune disorders in patients with vitiligo is approximately 10 times higher than the general population.¹⁰ The comorbidity including diabetes mellitus, thyroid disease, alopecia areata, RA, IBD, SLE and psoriasis.¹¹ They may share the same genetic susceptibility loci.¹² The types of the comorbidity may be influenced by race, area, sex, age and so on.^13,14 Psoriasis is the second most common comorbidity, especially in Asian patients.¹³

Psoriasis and vitiligo are two common autoimmune skin diseases that are prone to complications.¹⁵ On one hand, patients with vitiligo had an increased risk of psoriasis (fully adjusted HR, 1.71, 95% CI 1.35–2.17, p< 0.001). On the other hand, the incidence rate of new-onset psoriasis was estimated at 7.9 (95% CI 6.4–9.7) and 4.7 (95% CI 4.1–5.3) cases per 10,000 person-years among patients with vitiligo and control respectively.¹⁶ Patients with both vitiligo and psoriasis were older at the onset of the disease and had a greater prevalence of metabolic and cardiovascular comorbidities than those only with vitiligo.¹¹ In the pathogenesis of comorbidity, Th1 cytokines play an important role.¹⁷ They are both genetically susceptible with a high rate of family history and share common genetic susceptibility loci, such as rs9468925 in HLA-C/HLA-B.¹⁸ In recent years, immunotherapy has become an emerging therapy in various diseases. The immune checkpoint inhibitors have been used to treat multiple types of tumor,¹⁹ which trigger cutaneous immune-related adverse events however.²⁰ New onset or pre-existing exacerbation of psoriasis phenomenon were reported to be induced in oncology patients treated with PD-1 therapy,²¹ which may be associated with hyperactivation of effector T cells.²² At the same time, immunotherapy for melanoma often leaded to vitiligo-like depigmentation.²³ Patients with vitiligo have significantly shown more expression of PD-1 on CD8+ T cells compared with controls.²⁴

To explore the potential mechanisms underlying the comorbidity of two common chronic autoimmune skin diseases, we performed a differential analysis of two representative datasets of psoriasis (GSE30999)^25,26 and vitiligo (GSE75819)²⁷ stored in the GEO database. The DEGs obtained differentially expressed genes from GEO2R were subjected to functional enrichment analysis and PPI networks were constructed. Key modules and hub genes were extracted from the PPI network and functional enrichment analysis was performed. Then, psoriasis (GSE13355)^28,29 and vitiligo (GSE65127)³⁰ databases were selected from the GEO database for basic validation, and immunohistochemistry was performed at clinical samples for further validation. Finally, we paid attention to the immune checkpoint.

Materials and Methods

Data Mining and Analyze of GEO Database

The GSE30999 (psoriasis) and GSE75819 (vitiligo) datasets were obtained from the GEO database (https://www.ncbi.nlm.nih.gov/geo/). GSE30999 contains 170 expression datasets of moderate-to-severe psoriasis samples, including 85 samples of lesional skin and 85 samples of non-lesional skin. GSE75819 contains 30 expression datasets of vitiligo samples, including 15 samples of lesional skin and 15 of non-lesional skin. The online GEO2R tool was used to compare gene expression profiles between the diseased and control groups to determine the DEGs that must adjust both P-value< 0.05 and |log FC| (fold change) | ≥ 1. Probe sets without the corresponding gene symbols or genes with more than one probe set were excluded. Finally, common DEGs were identified using an online Venn diagram tool.

PPI Network Construction and the Selection of Hub Genes

The PPI of common DEGS was analyzed using the Search Tool for the Retrieval of Interacting Genes (STRING; https://cn.string-db.org/) (version 11.5). Interactions with a combined score over 0.4 were considered statistically significant and were then analyzed using Cytoscape (version 3.9.1). Single genes and gene pairs separated from the main body of the PPI network were removed. The online Metascape tool (http://metascape.org/) was used to analyze the key functional modules. Only physical interactions in STRING (physical score > 0.132) and BioGrid were used. Hub genes were identified using another plug-in, cytoHubba. Six common algorithms (MCC, MNC, Degree, EPC, Closeness, and Radiality) were used to evaluate and select the top 20 hub genes.

Gene Function Enrichment Analysis of the DEGS

Based on the Gene Ontology and biochemical pathway KEGG databases, the DEGs and hub genes were functionally enriched by cluster Profiler using Fisher’s exact test in R language. The online GeneMANIA tool (http://genemania.org/) was used to construct a co-expression network of hub genes to identify internal associations in gene sets.

Validation of Hub Genes Expression in Other Data Sets

The mRNA expression of the identified hub genes was verified using the GSE13355 and GSE65127 datasets obtained from the GEO database. The GSE13355 dataset contained 58 lesional skin samples from patients with psoriasis (PP), 58 non-lesional skin samples from patients with psoriasis (PN), and 64 normal skin samples from healthy controls (NN). GSE65127 consists of 10 lesional, 10 perilesional, and 10 non-lesional skin samples from patients with vitiligo and 10 control skin samples from healthy volunteers, which were referred to as VV, VP, VN, and NN, respectively. The Kruskal–Wallis rank-sum test was used for the statistical analysis between the two groups. Statistical significance was set at p < 0.05.

Immunohistochemical Staining Verification

Extracted skin tissues were embedded in OCT overnight and 5-μm -thick sections were obtained from frozen specimens. Tissue sections were fixed with 4% formaldehyde buffer and conducted antigen retrieval. The samples were quenched in 3% H₂O₂ for 10 min for endogenous peroxidase activity and treated with immunol staining blocking buffer (Beyotime, Shanghai, China) for 1 hour. Then the samples were incubated overnight at 4°C with primary antibodies against BUB1 (Proteintech, Wuhan, China, dilution 1:100), CDK1 (Proteintech, Wuhan, China, dilution 1:250). The secondary antibody (ZSGB-BIO, Beijing, China) for 30 minutes at 37°C and a DAB kit were applied to the samples. Micrographs of the stained sections were captured by light microscopy (Zeiss Imager A2, Germany), and three fields were randomly selected for each skin tissue section. The staining was scored according to the staining intensity and the area of the stained region. Intensity was evaluated as none (0), faint (1), moderate (2), strong (3). The sections were reviewed by two pathologists. The product of staining intensity and staining area scores was used as the H-score, which was evaluated using Mann–Whitney test. Statistical significance was set at p < 0.05.

Results

Identification of DEGs

Psoriasis and vitiligo are autoimmune diseases that share common drivers of pathogenesis. The comorbidity of these two diseases can sometimes been observed clinically, and the actual lesions are shown in Figure 1. Thus, we identified potential regulatory targets for both diseases using bioinformatic analysis. The flowchart of the study is presented in Figure 2. 3141 DEGs of psoriasis and 1100 DEGs of vitiligo were obtained with a threshold of |log FC| | ≥ 1 and P-value < 0.05 after standardizing the microarray results. The expression of DEGs is shown in the form of a volcanic map (Figure 3A and B). A total of 219 genes were obtained after the intersection of the psoriasis and vitiligo datasets (Figure 3C and 164 DEGs were obtained as they had the same expression trend (137 upregulated and 27 downregulated) (Figure 3D).

Figure 1 He skin lesion of the patient with psoriasis and vitiligo on left leg. (A). The right side of tibialis anterior. (B). The left side of tibialis anterior.

Figure 2 Research design flow chart.

Figure 3 Volcano diagram, Veen map and pie chart. (A). The volcano map of GSE30999. (B). The volcano map of GSE75819. (C). The Veen map of two datasets showed an overlap of 219 common DEGs. (D). The pie chart showed the trend constitution of common DEGs (137 up-regulated, 27 down-regulated, 55 contradictory and excluded genes.).

Analysis of Functional Characteristics of Common Genes

To explore probable pathways, common genes that focused on the GO and KEGG pathways were enriched. The enriched GO Biology Progress (BC) pathways contain “kinetochore organization”, “nuclear division”, “organelle fission” et al; the enriched GO Cellular Component (CC) pathways contain “chromosome, centromeric region”, “kinetochore”, “condensed chromosome, centromeric region” et al the enriched GO Molecular Function (MF) pathways contain “microtubule binding”, “tubulin binding”, “cysteine−type peptidase activity” (Figure 4A). In terms of KEGG Pathway, the four significant enrichment pathways are “NOD-like receptor signaling pathway”, “cell cycle”, “Epstein−Barr virus infection”, “RIG-I-like receptor signaling pathway” (Figure 4B).

Figure 4 Enrichment analysis of the 164 common DEGs with the same trend (A). GO enrichment analysis. (B). KEGG enrichment analysis. The size of the bar represents the number of genes involved, and the abscissa represents the frequency of the genes involved in the term total genes. The color of the bar represents p-value.

PPI Network Construction and Module Analyze

To further explore the interactions between common genes, the PPI network of the common DEGs with combined scores greater than 0.4 was constructed using Cytoscape (Figure 5). Genes were ranked in the form of concentric circles, according to the degree algorithm in the MCODE plugin. A larger size indicated a stronger relationship with other common genes. To explore the aggregation of genes and their corresponding functions, three main function modules were analyzed (Figures 6A–C) and the corresponding functions of each module are shown in Table 1. The main functions of the top 3 modules are “regulation of chromosome organization and separation”, “Interferon Signaling”, “nuclear chromosome segregation” respectively.

Table 1 The Matched Enriched Functions for the Top 3 Key Modules Analyzed from PPI Network Diagram of DEGs

Figure 5 PPI network diagram of 164 common DEGs with the same expression trend.

Figure 6 The top 3 key functional modules analyzed from PPI network diagram of DEGs. (A). The gene clustering module “MCODE_1” which is correlated with the function “regulation of chromosome organization and separation”. (B). The gene clustering module “MCODE_2” which is correlated with the function “interferon signaling”. (C). The gene clustering module “MCODE_3” which is correlated with the function “nuclear chromosome segregation”.

Selection and Analyze of Hub Genes

Using the six algorithms of the cytoHubba plugin in Cytoscape, we identified the top 20 hub genes (Table 2). After taking the intersection of the Venn diagrams, we found 14 common hub genes (Figure 7A) (BUB1, CEP55, CDK1, TOP2A, CENPF, PBK, MELK, CCNB2, MAD2L1, NUSAP1, TTK, NEK2, CDKN3, PTTG1) and Table 3 shows their full names and related functions by GeneCards.

Table 2 The Top 20 Hub Genes Ranked in cytoHubba

Table 3 The Details of the 14 Overlapping Hub Genes

Figure 7 Selection and Analyze of Hub Genes. (A). Veen map. The six algorithms showed an overlap of 14 hub genes. (B). Co-expression network of hub genes analyzed via GeneMANIA. (C).GO Enrichment analysis of the hub genes. D.KEGG enrichment analysis of the hub genes. The size of the bar represents the number of genes involved, and the abscissa represents the frequency of the genes involved in the term total genes. The color of the bar represents p-value.

Based on the GeneMANIA database, we analyzed the co-expression network and related functions of the hub genes. These genes showed a complex PPI network with 78.91% co-expression, 8.74% physical interactions, 5.83% pathways, 4.04% colocalization, 2.46% predicted, and 0.02% Genetic Interactions (Figure 7B).

The enriched GO BC pathways contain “nuclear division”, “organelle fission”, “sister chromatid segregation” et al; the enriched GO CC pathways contain “chromosomal region”, “chromosome, centromeric region”, “condensed chromosome” et al; the enriched GO MF pathways contain “protein serine kinase activity”, “protein serine/threonine kinase activity”, “transmembrane receptor protein tyrosine kinase activity” et al (Figure 7C). In terms of KEGG Pathway, the three significant enrichment pathways are “cell cycle”, “oocyte meiosis”, “progesterone-mediated oocyte maturation” et al (Figure 7D).

Validation of Hub Genes Expression in Other Data Sets and at the Clinical Sample

To verify the reliability of these hub gene expression levels, we selected two other data series (GSE13355 and GSE65127) containing psoriasis and vitiligo lesions and analyzed the expression levels of these hub genes. The results showed that compared with normal skin and non-lesional skin, all 14 hub genes were significantly upregulated in psoriatic skin lesions (Figure 8). However, only seven hub genes were validated as highly expressed in the vitiligo dataset (Figure 9), including BUB1, CDK1, PBK, CCNB2, MAD2L1, TTK, and NEK2, indicating that these seven genes should be given more attention. Validation of hub genes in the vitiligo data series was limited by the paucity of publicly available vitiligo datasets for analysis and small number of samples in the available vitiligo data series. Therefore, the remaining seven genes should also be considered, including CEP55, TOP2A, CENPF, MELK, NUSAP1, CDKN3, and PTTG1.

Figure 8 The expression level of hub gene in GSE13355. The comparison between two datasets used the mean Kruskal–Wallis rank sum test. P-value < 0.05 was considered statistically significant. PP, lesional skin samples from psoriatic patients; PN, non-lesional skin samples from psoriatic patients; NN, normal skin samples from healthy controls. * p < 0.05; ** p < 0.01; ***p < 0.001; ****p < 0.0001.

Figure 9 The expression level of hub gene in GSE65127. The comparison between two datasets used the mean Kruskal–Wallis rank sum test. P-value < 0.05 was considered statistically significant. LL, lesional skin samples from vitiligo patients; LN, non-lesional skin samples from vitiligo patients; NN, control skin samples from healthy volunteers. * p < 0.05; ** p < 0.01; ***p < 0.001.

BUB1 is a kind of serine/threonine-protein kinase, which is essential for spindle-assembly checkpoint signaling and for correct chromosome alignment during mitosis. It has a key role in the assembly of checkpoint proteins at the kinetochore, being required for the subsequent localization of CENPF, BUB1B, CENPE and MAD2L1. CDK1 plays a key role in the control of the eukaryotic cell cycle by modulating the centrosome cycle as well as mitotic onset. To further validate the feasibility and authenticity of the screened hub genes, six comorbidity samples and six healthy samples were collected from the patients of our department and performed IHC staining. The semiquantitative H-score of BUB1, CDK1 in comorbidity of psoriasis and vitiligo skin tissue were higher than those in HC skin tissue (p< 0.05) (Figure 10A–C). It is agreed with the findings of psoriasis (GSE30999) and vitiligo (GSE75819), and the validation of psoriasis (GSE13355) and vitiligo (GSE65127).

Figure 10 The expression level of hub gene at clinical sample. (A). Immunohistochemical staining for BUB1 and CDK1 in the lesional skin samples from patients with the comorbidity of psoriasis and vitiligo, and in the control skin samples from healthy volunteers (10×, 40× magnification). Scale bars = 200 and 50 μm for the low (10×) and high magnification pictures (40×), respectively. (B). The H-score of BUB1 for two groups presented by boxplot with Mann–Whitney test. (C). The H-score of CDK1 for two groups presented by boxplot with Mann–Whitney test. P-value < 0.05 was considered statistically significant. CO: lesional skin samples from patients with the comorbidity of psoriasis and vitiligo; HC, control skin samples from healthy volunteers. ** p < 0.01.

The Common Immune Checkpoint Between Psoriasis and Vitiligo

Immune checkpoints is always considered to play an important role in the development of cancer and other chronic infectious diseases. However, in psoriasis, many checkpoints are differentially expressed, including CD274, CTLA4, HAVCR2, LAG3, PDCD1, PDCDILG2, SIGLEC15, and TIGIT. CTLA4 is highly expressed in vitiligo. So, Cytotoxic T lymphocyte-associated antigen-4 (CTLA4) may be the probable common immune checkpoint in disease pathogenesis of comorbidity (Figure 11).

Figure 11 The Common Immune Checkpoint between Psoriasis and Vitiligo. (A). The differentially expressed immune checkpoint in psoriasis. (B). The differentially expressed immune checkpoint in vitiligo. ** p < 0.01; ***p < 0.001.

Discussion

Bioinformatics analysis of the disease is currently a common investigation method. Previous studies have separately explored hub genes associated with psoriasis or vitiligo. However, few studies have focused on analyzing the analyze uncovering of co-morbidities of psoriasis and vitiligo using advanced bioinformatic methods. Owing to the interesting and amazing clinical characteristics of these two diseases, we identified the combined genetic profiles, including common DEGs and hub genes, and evaluated the potential mechanism of the occurrence of this complication explored for the first time, which helped to further clarify the potential mechanism of the comorbidity.

We identified 219 overlapping DEGs between psoriasis and vitiligo, of which 164 showed the same expression trend (137 upregulated and 27 downregulated). Enrichment analysis of DEGs significantly focused on the cell cycle and mitotic cell division, involving kinetochore organization and nuclear fission, and on the role of cytoskeletal reorganization and chromatin remodeling, including the microtubule binding and tubulin binding pathways. Among 14 top screened hub genes, seven hub genes (BUB1, CDK1, PBK, CCNB2, MAD2L1, TTK, and NEK2) were validated high expression in other psoriasis and vitiligo datasets, but the remaining seven hub genes (CEP55, TOP2A, CENPF, MELK, NUSAP1, CDKN3, and PTTG1) were validated as being highly expressed only in other psoriasis datasets. It does not indicate that these genes are not significant, but is limited by the lack of publicly available data on the web. This is because vitiligo can be easily diagnosed by clinical symptoms, and so few patients can be receptive to pathologic sampling. Thus, small number of samples in the available vitiligo data series online can be used.

The hub genes were divided into the following three categories: genes associated with mitotic spindle formation and checkpoints, including MAD2L1, TTK, BUB1, NEK2, and CENPF; genes associated with cell cycle regulation, including PBK, MELK, CEP55, CCNB2, CDK1, and CDKN3; and genes associated with microtubule-binding proteins, including NUSAP1, DNA topoisomerase TOP2A, and PTTG1, which maintain chromosome stability.

Based on these results, hub genes showed a strong connection with spindle localization and the attachment of spindle microtubules to kinetochores. The cytoskeleton is the protein fiber reticulum system in eukaryotic cells and consists of tubulin-containing microtubules, actin-containing microfilaments, and intermediate filaments. It is associated with maintenance of cell shape and tissue integrity. It can resist external stress, modulate cellular motility, regulate intracellular microfilament movement, and facilitate the onset of signal transduction.³¹ Traditional studies on immune diseases have been limited to these pathways. However, the effect of concentration on the mechanical properties of the cells remains unclear. Mechanical force not only occurs at the level of macroscopic objects but also in cells and molecules.³² Previous studies have shown that chimeric antigen receptor (CAR) T cells function by coordinating multiple mechanical forces during immune disorders. CAR T-cell activation is a mechanotransductive process that involves spatiotemporal remodeling of the cytoskeleton.^33,34 It may affect the ion channel, metabolic efficiency, and genetic characteristics of CAR cells, and modulate immune function.³²

The eukaryotic cell cycle consists of the G1, S, G2, and M phases. In human cell culture, the interphase of mitosis occupies 23 of its 24 h, including the G1, S, and G2 phases. Ordinary cyclin proteins are expressed during interphase. CDKs, which vary cyclically, initiate the cell cycle machinery.³⁵ Among these, CDK1 strictly controls cell cycle phases and mitotic events. CCNB2 is a monitor protein that plays an important role in the switch from G2 to M phase. A previous study has shown that it may promote the release of key molecular targets of psoriasis by regulating mast cell activation and macrophage polarization.³⁶ TTK is another important cell cycle regulator. BUB1 is a part of the mitotic checkpoint complex (MCPC).³⁷ Checkpoints are biological processes that allow the cell cycle to halt and can be used to monitor DNA stability. The MCPC is fundamental to the spindle assembly checkpoint (SACP). A previous study has demonstrated the important role of BUB1 in psoriasis.³⁸ A lack of BUB1B can result in aneuploidy and chromosomal instability, which may increase the risk of cancer. Genetic alterations in BUB1B in psoriasis are closely associated with immune cell infiltration.³⁹ As an SACP, MAD2L1 ensures that chromosomes are appropriately oriented toward the metaphase plate during cell division.⁴⁰ CENPF is a p53 signaling pathway protein whose synthesis has a positive correlation with the number of CD8⁺ T-cells and a negative correlation with the number of CD4⁺ T-cells.^41,42 CTLA-4 is an immune checkpoint protein that is transported to the cell surface upon T-cell activation from intracellular recycling vesicles⁴³ and is responsible for the transmission of negative regulation signals. The imbalance of CTLA-4 leads to the overactivation of CD8+ T-cells, which in turn attack the melanocytes and cause vitiligo.⁴⁴ Abatacept is a selective T-cell co-stimulatory modulator, a fusion protein consisting of the IgG Fc segment and the extracellular structural domain of CTLA-4, which can inhibit T-cell activation. Clinical trials of Abatacept for the treatment of vitiligo are currently underway.⁴⁵ The CTLA-4 signaling peptide showed a significant reduction of psoriatic skin inflammation with increased Treg cell proportion and reduced IL-17 production by T cells, indicating a potential role in modulating psoriatic skin disease.^1,46 Thus, the imbalance of CTLA-4 may play a key role in the immune treatment of comorbidity in dermatology.

Based on these findings, regulation of cytoskeletal dysfunction involved in mitotic cell division may play an important role in the pathogenesis of comorbidity of psoriasis and vitiligo. It may function in both spatial and temporal dimensions. In the spatial dimension, the cytoskeleton of cells is an important factor in nuclear division, organelle fission, and sister chromatid separation, with dynamic involvement of microtubules and microfilaments in the cell machinery. In the temporal dimension, the checkpoints and phases of the cell cycle are significant for mitotic entry, mitosis, or mitotic exit.

Conclusions

In summary, we selected seven hub genes, BUB1, CDK1, PBK, CCNB2, MAD2L1, TTK, and NEK2, that may be involved in the pathogenesis of both psoriasis and vitiligo. The mitotic cell division, cell cycle and the cytoskeleton could play crucial roles in the pathogenesis of comorbidity. And we proposed a new spatio-temporal theory for the pathogenesis of the comorbidity, cytomechanics-guided cytoskeleton remodeling as the spatial dimension, and the cell cycle and mitosis as the temporal dimension. The disruption of these two dimensions may play an important role in psoriasis-vitiligo comorbidities pathogenesis, which has been overlooked in the previous studies. It is also proposed for the first time that highlight the function of immune checkpoint in the pathogenesis of the comorbidity of these two diseases. CTLA-4 may become a new treatment target for the further research.