Identification of key genes and evaluation of immune cell infiltration in vitiligo

: Background: To improve the understanding of the molecular mechanism of vitiligo is necessary to predict and formulate new targeted gene therapy strategies. Methods: GSE65127, GSE75819, GSE53146 and GSE90880 were collected, and obtained four groups of differentially expressed genes (DEGs) by limma R package. Through weighted gene co-expression network analysis (WGCNA), the co-expression of genes with large variance in GSE65127 and GSE75819 was identified. Enrichment analysis of intersection gene between module genes and DEGs with the same up-regulated or down-regulated in GSE65127 and GSE75819 was performed. In addition, ssGSEA was used to identify the immune infiltration of vitiligo in four datasets. Results: A total of 3083 DEGs and 16 modules were identified from GSE65127, and 5014 DEGs and 6 modules were screened from GSE75819. Finally, 77 important DEGs were identified. Enrichment analysis showed that 77 DEGs were mainly involved in spliceosome etc. The results of GSVA showed that melanogenesis, Fc gamma R-mediated phagocytosis, leishmaniasis, Wnt pathway and glycolipid metabolism were important KEGG pathways. The genes involved in these pathways were identified as key genes (MARCKSL1, MC1R, PNPLA2 and PRICKLE2). The AUC values of MC1R were the highest. Furthermore, different immune cells had different infiltration in vitiligo. There was a high correlation between immune cells and key genes. Conclusions: MC1R was found as a key gene in vitiligo and involved in the melanogenesis. The immune cells were different infiltration in vitiligo. These results suggested that key genes may be used as markers of vitiligo, and were associated with immune cell, especially MC1R.


Introduction
Vitiligo is a multifactorial and complex skin disease characterized by flaky epidermal skin peeling [1]. The prevalence rate in the world population is 0.5-1%, and there is no obvious gender bias [2]. The highest prevalence was reported in India (8.8%), followed by Mexico (2.6-4%) and Japan (1.68%) [3]. Because of their appearance and subjective symptoms, patients with dermatosis often suffer from psychological stress and quality of life [4]. Research showed that patients with vitiligo have emotional disorders such as embarrassment, anger, worry or frustration [5]. The incidence rate of vitiligo is often manifested by obvious family clustering. It has been reported that 20% of vitiligo patients also have relatives diagnosed with vitiligo [6].
It has been recognized that predisposing factors, including exposure to sunlight or skin trauma, lead to melanocyte oxidative stress and T cell-mediated autoimmune response [7]. In fact, the etiology and pathogenesis of vitiligo are not fully understood. However, the loss and dysfunction of melanocytes are the main causes of vitiligo [8]. Melanocyte disappearance may involve a variety of mechanisms, including genetic susceptibility, environmental incentives, impaired melanocyte regeneration, and changes in inflammatory response [9]. Autoimmune and oxidative stress are also two important mechanisms of vitiligo [10]. In addition, metabolic abnormalities such as insulin resistance, dyslipidemia and skin involvement can be observed in vitiligo [11].
However, there are currently no FDA-approved medical treatments for vitiligo, and available off-label treatments are problematic and often ineffective [12,13]. The existing treatment methods, such as phototherapy, local immunomodulators and surgical methods, partially solve the problem of immune and melanin reduction [14]. There are also quite a number of patients who do not get enough help, or the success of treatment only lasts for a short time [15]. Therefore, the development of safe and effective therapies requires a better understanding of the pathogenesis of the disease to identify new therapeutic targets.
This study was based on bioinformatics method to identify the key target genes of vitiligo patients. Furthermore, the potential mechanism of target gene was explored by recognizing the immune infiltration. This will further deepen our understanding of the mechanism of vitiligo, and provide more and wider options for clinical treatment.

Data collection
We collected mRNA expression profiles of Vitiligo from the Gene Expression Omnibus (GEO) database. Expression profiling of skin tissue by array in GSE65127 included 10 vitiligo patients and 10 healthy volunteers from France. There were 30 samples from 15 individuals' lesional and non-lesional skin of vitiligo patients came from India in GSE75819. GSE53146 included skin samples of 5 vitiligo patients and 5 from controls from USA. GSE90880 included whole blood samples of 8 vitiligo patients and 6 healthy controls from USA. The data were normalized with the normalizebetweenarrays function of the limma package.

Identification of differentially expressed genes
The differentially expressed genes (DEGs) between vitiligo and control were obtained with limma R package [16,17]. P < 0.05 was considered as significant difference [18].

Construction of co-expression network
The co-expression network of first quarter of genes with large variance was constructed by weighted gene co-expression network analysis (WGCNA). The optimal soft threshold of adjacency calculation is determined by graphic method. The transformed expression matrix is input into WGCNA R package [19] functions and modules, and the eigengenes are obtained.

Enrichment analysis
Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis was performed by clusterProfiler R package [20,21] to identify the biological function. In addition, the KEGG pathway involved by genes also was identified by gene set enrichment analysis (GSEA) software. P < 0.05 was set as the cut-off criterion.
The results of gene enrichment were quantified using gene set variation analysis (GSVA) R package. GSVA scores were calculated using a Kolmogorov-Smirnoff-like random walk statistic and a negative value for a particular sample and gene set. Using limma R package for differential analysis could find KEGG pathway with significant difference among samples.

Infiltration of immune cells
The marker gene set for immune cell types was obtained from Bindea et al. [22]. We used the single-sample gene set enrichment analysis (ssGSEA) program to derive the enrichment scores of each immune-related term. The infiltration levels of immune cells were quantified by ssGSEA R package. Pearson correlation was used to calculate the correlation among immune cells or between immune cells and genes.

Important differentially expressed genes in vitiligo
By comparing the differences between vitiligo patients and controls, 3083 differentially expressed genes (DEGs) were found in GSE65127 ( Figure 1A, Table S1). The co-expression analysis of genes with large variance showed that a total of 5014 genes had synergistic expression effect and were clustered into 16 co-expression modules ( Figure 1B). In addition, 8264 DEGs were obtained in GSE75819 ( Figure 1C, Table S2). A total of 4092 genes with large variance were selected and clustered into 6 co-expression modules ( Figure 1D). In order to identify the important DEGs related to vitiligo, we screened 77 DEGs from the two groups of co-expression module genes that were up-regulated and down-regulated simultaneously in two datasets ( Figure 1E). The clustering of 77 DEGs may distinguish a subset of vitiligo and control samples in GSE65127, but clustering in GSE75819 may completely distinguish vitiligo and control samples. Since the expression of these genes was significantly different between the two groups of samples ( Figure 2).

Biological function and signal pathway
The enrichment analysis of 77 DEGs showed that they were significantly involved in 268 biological processes (BP). It mainly included response to light stimulus, response to anticoagulant, AMP metabolic process, peroxisome localization and mitotic nuclear division ( Figure 3A). KEGG enrichment results showed that the genes were significantly involved in five terms, including spliceosome, Nucleotide excision repair, Cell cycle, Lysosome and Riboflavin metabolism ( Figure  3B). In the GSEA results of the two datasets, there were three identical terms, including spliceosome, RNA transport and basal transcription factors ( Figure 3C).

Key pathways and genes in vitiligo
The GSVA method was used to quantify the expression of KEGG pathway of 77 DEGs in two datasets ( Figure 4A). The up-regulated or down-regulated signaling pathways in both datasets were obtained. The top five largest |logFC| pathways were considered as closely related to vitiligo ( Figure  4B). MC1R, MARCKSL1, PNPLA2 and PRICKLE2 were involved in these five pathways ( Figure  4C). They were down regulated in vitiligo ( Figure 4D). ROC curve showed that their AUC values were greater than 0.8 in both GSE65127 and GSE75819 ( Figure 4E,F). In particular, the ROC values of MC1R were the highest.

Difference of immune infiltration in vitiligo
Comparing the difference of 24 kinds of immune cells between vitiligo and control in GSE65127, GSE75819, GSE53146 and GSE90880, it was suggested that immune cells were involved in the pathogenesis of vitiligo ( Figure 5A). Among them, Th2 cells and aDC were up-regulated in four datasets, while B cells were down-regulated. These immune cells were clustered into four groups, and there was a certain correlation between different types of cells ( Figure 5B). The positive correlation between cytotoxic cells and NK cd56dim cells was the strongest ( Figure 5C). Among these immune cells, NK cd56bright cells had the strongest positive correlation with MC1R, T cells and Tgd had the strongest negative correlation with PRICKLE2 ( Figure 5D).

Discussions
Treatment and management of vitiligo has always been a tough challenge for scientists and dermatologists. Therefore, vitiligo may be treated as a chronic immune disease [23]. In this study, we obtained DEGs related to vitiligo by analyzing the transcriptome data of vitiligo in the public database. The key pathways of inflammation and metabolism were identified in enrichment analysis. Moreover, the genes involved in these pathways were also identified as key target genes and have high correlation with immune cells.
Among the four key genes identified in this study, melanocortin 1 receptor (MC1R) is a G protein coupled receptor that regulates melanocyte proliferation and function [24]. MC1R signal enhanced melanin synthesis, induced antioxidant activity, enhances DNA repair process, and regulates inflammation [25]. MC1R was also considered to play a role in the pathogenesis of vitiligo [26]. Consistent with our analysis results, qPCR assay also confirmed that MC1R expression was significantly low in vitiligo patients [27]. Myristoylated alanine-rich C kinase substrate-like 1 (MARCKSL1) showed antiangiogenic effect by inhibiting the phosphorylation of VEGFR-2 dependent Akt / PDK-1 / mTOR [28]. In patients with vitiligo, angiogenesis and mast cells increased at the sites with low lymphocyte count, which may be one of the reasons for secondary inflammatory process of vitiligo [29]. The loss of patatin-like phospholipase domain containing 2 (PNPLA2) increased the cell resistance to oxidative stress [30]. Oxidative stress can induce the formation of autoimmunity in vitiligo [31]. It had been reported that mutations in the PRICKLE2 gene were associated with human diseases such as epilepsy and autism spectrum disorders [32]. Our results suggested that PRICKLE2 may inhibit the progression of vitiligo, but there is no other study on its relationship with vitiligo.
Among the biological functions and signal pathways we analyzed, melanogenesis has been considered as a key factor in the formation of vitiligo [14,33]. Various therapeutic methods for melanogenesis have been studied and discussed [34,35]. Fc gamma R (FcγR) -mediated phagocytosis has the ability to activate the immunoregulation pathway [36]. Wnt signal regulation may promote appropriate melanocyte regeneration in vitiligo [37]. Decreased Wnt activation played a role in preventing melanocyte differentiation in pigmented vitiligo [38].
More evidence showed that the dynamic and pathogenesis of vitiligo was closely related to the immune response of skin and melanocytes. Studies had shown that anti melanocyte specific cytotoxic T cells played a central role in the final effector phase [39]. In fact, vitiligo is usually associated with autoimmune diseases; genome wide association studies and functional pathway analysis show that most of vitiligo susceptible sites encode components of the immune system; immune cells are found at the edge of active peeling in patients with vitiligo [40].
As with all studies, our study on the molecular mechanism of vitiligo has limitations. First of all, our analysis data was from the public database, and the analysis results need to be verified by molecular experiments. Secondly, the molecular mechanism related to the key genes of vitiligo need to be further verified by a series of experiments. Finally, whether the selected key genes and molecular mechanisms can be used as therapeutic targets of vitiligo also need a large number of samples and clinical trials to verify.
In conclusion, by combining WGCNA with analysis of differentially expressed gene, our study had generated potential diagnostic value of treatment-related genes and molecular mechanisms. The new treatment methods may be based on the in-depth understanding of the pathogenesis of the disease to adopt more targeted methods, which may provide higher efficacy.