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

1.   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.

2.   Materials and methods

2.1.   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.

2.2.   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 ]} .

2.3.   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.

2.4.   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.

2.5.   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.

3.   Results

3.1.   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 ).

Figure 1.  Screening of important differentially expressed genes in vitiligo. (A) Volcano map of differentially expressed genes between vitiligo and control in GSE65127. (B) The co-expression modules of genes with large variance in GSE65127. (C) Volcano map of differentially expressed genes between vitiligo and control in GSE75819. (D) The co-expression modules of genes with large variance in GSE75819. (E) The intersection of co-expression module genes and differentially expressed genes that were up-regulated or down regulated simultaneously in two datasets.

DownLoad: Full-Size Img PowerPoint

Figure 2.  Cluster heatmap of important differentially expressed genes in two sets of data.

DownLoad: Full-Size Img PowerPoint

3.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 ).

Figure 3.  Important DEGs are involved in biological functions and signaling pathways. (A) The biological processes of important DEGs. (B) The KEGG pathway of important DEGs. (C) The same KEGG pathway in GSEA of two datasets for important DEGs.

DownLoad: Full-Size Img PowerPoint

3.3.   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.

Figure 4.  Key signaling pathways and genes related to vitiligo. (A) The signal pathways were up-regulated or down regulated simultaneously in two datasets. (B) The orders of | logFC | for signal pathways. (C) A network of genes involved in signaling pathways. (D) Genes involved in the top five signaling pathways. (E) ROC curve of four genes in two sets of data. (F) The expression value and AUC value of four genes in two sets of data.

DownLoad: Full-Size Img PowerPoint

3.4.   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 ).

Figure 5.  Immune status of vitiligo. (A) The difference of immune infiltration between vitiligo and control. (B) The immune cells were clustered into four different categories. (C) Correlation between immune cells in vitiligo or control. (D) Correlation between immune cells and key genes.

DownLoad: Full-Size Img PowerPoint

4.   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.

Acknowledgments

This work was funded by Tianshan Cedar Project (2018XS20).

Data availability statement

The gene expression files were obtained from GEO database: GSE65127, GSE75819, GSE53146 and GSE90880.

Conflict of interest

The authors declare no conflict of interest in this paper.