Identification of key genes and pathways involved in vitiligo development based on integrated analysis.

Vitiligo is a chronic skin condition lack of melanocytes. However, researches on the aetiology and pathogenesis of vitiligo are still under debate. This study aimed to explore the key genes and pathways associated with occurrence and development of vitiligo.Weighted gene coexpression network analysis (WGCNA) was applied to reanalyze the gene expression dataset GSE65127 systematically. Functional enrichments of these modules were carried out at gene ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), gene set variation analysis (GSVA), and gene set enrichment analysis (GSEA). Then, a map of regulatory network was delineated according to pivot analysis and drug prediction. In addition, hub genes and crucial pathways were validated by an independent dataset GSE75819. The expressions of hub genes in modules were also tested by quantitative real-time polymerase chain reaction (qRT-PCR).Eight coexpressed modules were identified by WGCNA based on 5794 differentially expressed genes of vitiligo. Three modules were found to be significantly correlated with Lesional, Peri-Lesional, and Non-Lesional, respectively. The persistent maladjusted genes included 269 upregulated genes and 82 downregulated genes. The enrichments showed module genes were implicated in immune response, p53 signaling pathway, etc. According to GSEA and GSVA, dysregulated pathways were activated incessantly from Non-Lesional to Peri-Lesional and then to Lesional, 4 of which were verified by an independent dataset GSE75819. Finally, 42 transcription factors and 228 drugs were spotted. Focusing on the persistent maladjusted genes, a map of regulatory network was delineated. Hub genes (CACTIN, DCTN1, GPR143, HADH, MRPL47, NKTR, NUF2) and transcription factors (ITGAV, SYK, PDPK1) were validated by an independent dataset GSE75819. In addition, hub genes (CACTIN, DCTN1, GPR143, MRPL47, NKTR) were also confirmed by qRT-PCR.The present study, at least, might provide an integrated and in-depth insight for exploring the underlying mechanism of vitiligo and predicting potential diagnostic biomarkers and therapeutic targets.

Introduction

Vitiligo is a common acquired, chronic skin condition with the loss of functioning melanocytes, affecting skin, hair or mucosa, or all, characterized by milky white, nonscaly spots with distinct margins.[ These clinical symptoms seriously affect the mental health and quality of life of individuals.[ With an estimated prevalence of 0.5% to 1% in the world population,[ vitiligo is regarded as one of the most important skin conditions causing prejudice and segregation in certain cultures. Both sexes in adults and children seem to be affected by vitiligo equally.[ However, the etiology and pathogenesis of vitiligo is still obscure and both multifactorial and polygenic. Selective destruction of the pigment-producing cell, namely melanocytes, is core to vitiligo.[ Autoantibodies against melanocytes were shown to exist in patient serum, and lately have been considered as predictors of disease progression.[ In addition, several autoimmune comorbidities have been reported in association with vitiligo such as autoimmune thyroid disease,[ systemic lupus erythematosus,[ Addison disease[ and psoriasis,[ etc. Furthermore, a number of susceptibility loci for vitiligo were detected through genome-wide association studies, lots of which are shared with other autoimmune diseases.[

To regain and maintain the pigment, current therapeutic options for vitiligo are mainly based on 3 aspects: modulating the autoimmune response; reducing oxidative stress in melanocytes; activating regeneration of melanocytes.[ However, affected individuals require long-term treatment and poor efficacy. So far, vitiligo remains a disorder with restricted therapeutic options and difficult to control its progression effectively for patients, which underlines the unmet need to define the molecular mechanism underlying this disorder and to clear the therapeutic target for new safe and effective therapies to manage this skin condition. Therefore, more intensive research is needed to discover the molecular mechanism, detect novel promising therapeutic targets, and develop efficient therapeutic approaches for vitiligo.

Currently, insights into the pathogenesis of autoimmune diseases at the system level can be provided by researches of systems biology, which can facilitate biological studies. Some datasets of genome-wide gene expression in patients with vitiligo are freely accessible in the public databases, viz. Gene Expression Omnibus.[ Much valuable information in those datasets can be mined and reused by some tools of systems biology. Consequently, weighted gene coexpression network analysis (WGCNA) was performed to categorize the genes expressed aberrantly and significantly into several functional modules,[ expecting to deepen insights into the pathogenesis of vitiligo and seek some novel promising drug targets for the therapy of vitiligo.

Materials and methods

Data resources and principal component analysis (PCA)

The vitiligo datasets GSE65127 and GSE75819 were collected from the NCBI Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/).[ GSE65127 includes a total of 40 samples obtained from 10 patients and 10 healthy volunteers (1 biopsy in matched anatomical areas), comprising tissues from 10 Lesional skins, 10 Peri-Lesional skins, 10 Non-Lesional skins, and 10 healthy volunteers’ skins. GSE75819 contains tissues from 15 Lesional skins and 15 Non-Lesional skins. Series matrix files and clinical information tables were downloaded from the GEO website. PCA is an unsupervised method and data reduction technology, which allows the analysis of major sources of change in multidimensional datasets without introducing inherent deviations. Hence, PCA were performed on disease and health samples to condense the information down to a reduced number of linearly independent variables.

Differential expression analysis

The “sva” R package was applied to conduct batch normalization of the expression data from the 3 different datasets. Then, the normalized expression matrixes of Non-Lesional, Peri-Lesional, and Lesional differences were analyzed by “limma” R package respectively.[ The differentially expressed genes (DEGs) were acquired by setting threshold P < .05.

Coexpression analysis

WGCNA is a method to establish scale-free gene coexpression network.[ A weighted gene coexpression network of 3 groups of DEGs was constructed using “WGCNA” R package. The soft threshold power of β was set as 9, and weighted adjacency matrix was subsequently generated. Then, hierarchical clustering analysis was carried out. And a weighted adjacency matrix was formed. In addition, the weighted adjacency matrix was transformed into a topological overlap matrix to estimate its connectivity in the network. To evaluate the association between gene coexpression modules and traits, “nlme” R package was adopted to establish a linear mixed model. For each WGCNA module, the gene with kme Pearson correlation (> 0.90)[ was regarded as a hub gene. Afterward, the receiver operating characteristic (ROC) curves of hub genes were analyzed by “pROC” R package.[

GO function and KEGG pathway enrichment analysis and gene set enrichment analysis (GSEA)

To clarify the possible biological roles of these genes in coexpression networks, the cluster Profiler R package[ was performed to create gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment paths, results, and plots. P < .05 indicated GO and KEGG pathway with significant difference. Thereafter, GSEA[ was employed according to the application available at the Broad Institute Gene Set Enrichment Analysis website. For purpose of testing the KEGG pathways gene sets, the activation or inhibition of KEGG pathway was also evaluated by GSVA algorithm.[

Pivot analysis and drug prediction

Pivot is defined as a regulator, including transcription factors, which has a significant regulatory effect on module genes. To explore the pivotal regulators of vitiligo-related module genes, all human transcription factor target data in TRRUST V2 database[ was utilized as background set for pivot analysis. Setting threshold P < .05. Prediction of drug and drug targets for the vitiligo dysfunction module were also carried out. On the basis of drug target data from Drugbank database,[ drugs with significant regulatory modules were screened according to hypergeometric test. Setting threshold P < .05.

RNA expression analysis by quantitative reverse transcription polymerase chain reaction (qRT-PCR)

The whole blood samples of 3 vitiligo patients and 3 healthy controls were collected respectively. Written informed consent was obtained from each subject, and the research protocol was approved by the ethics committee of People's Hospital of Xinjiang Uygur Autonomous Region. The total RNA of whole blood samples were extracted using Trizol reagent (Invitrogen, USA) by following the manufacturer's protocol. Following DNase treatment, RNA was reverse transcribed into complementary DNA with a complementary DNA synthesis kit (Takara, Dalian, China), according to the manufacturer's protocol. qPCR was performed by using SYBR Select Master Mix (Applied Biosystems, USA) and primers (Table 1). The expression of target gene was normalized to GAPDH levels in respective samples as an internal control and calculated using the 2−ΔΔCT method. Students t test was used to compare the differences between 2 groups. P < .05 was considered statistically significant.

Table 1

Primer sequence.

Results

Coexpression of differentially expressed genes in vitiligo

First, PCA was used to analyze the degree of dispersion of different samples. It was found that the 3 groups of disease samples were closer, while the healthy control group samples were farther (Supplemental Digital Content (Figure S1)). Then, DEGs were screened between Non-Lesional, Peri-Lesional, Lesional, and healthy skin tissues, respectively, which were considered to be associated with vitiligo. 3132 DEGs were identified in Non-Lesional, 3610 in Peri-Lesional, and 3083 in Lesional (Fig. 1A). Next, WGCNA was performed on the 3 groups of DEGs and a gene coexpression network was constructed. By eliminating the gray module without cooperative expression behavior, 8 modules were obtained, in which the genes had the behavior of coexpression (Fig. 1B–D). Through the correlation study between module and phenotype, it was noteworthy that the highest correlation with Lesional, Peri-Lesional, and Non-Lesional was module brown, blue, and green respectively (Fig. 1E). In addition, gene with the highest connectivity, named hub gene, were identified in each module (Table 2). Hub gene plays a significant role in the module and may function as a key part in the disease. As well, by ROC curve, the ability of hub genes to determine the progression of vitiligo was evaluated (Fig. 2A–H). The Area Under Curve value of NUF2 and HSPB6 increased gradually from Non-Lesional to Lesional, indicating that the more severe the degree of vitiligo, the stronger the evaluation ability of NUF2 and HSPB6. Surprisingly, the diagnostic ability of 7 hub genes was verified in an independent gene expression dataset GSE75819 (Fig. 2I, J). These results implied that most of the DEGs were coexpressed and interacted with each other in the occurrence and development of vitiligo.

Figure 1

Coexpression of differentially expressed genes in vitiligo was observed. A, The differentially expressed genes between the skin tissue of patients with vitiligo in Lesional, Peri-Lesional, Non-Lesional, and the healthy controls. B, Soft threshold power analysis was used to obtain the scale-free fit index of network topology. C, Hierarchical cluster analysis was conducted to detect coexpression clusters with corresponding color assignments. Different colors represent different modules. D, Heatmap depicts the Topological Overlap Matrix (TOM) of genes selected for weighted coexpression network analysis. Light color represents lower overlap and red represents higher overlap. E, The correlation between module and phenotype. Red reflects positive correlation, blue reflects negative correlation.

Table 2

Hub genes in modules of differentially expressed genes in vitiligo.

Figure 2

ROC curves of hub genes. A–H, Calculation of Area Under Curve value of hub genes in the modules and ROC curve drawing. Different colors represent different groups. I, J, The seven hub genes validated in dataset GSE75819 and corresponding Area Under Curve values. ROC = receiver operating characteristic.

Molecular dysregulation mechanism in vitiligo

These module genes may play different roles in regulating the occurrence and development of vitiligo. To identify the biological functions of dysfunctional genes, GO and KEGG enrichment analysis of module genes were performed. Abundant GO terms were collected, including 562 cell composition terms, 842 molecular functional terms, as well as 3916 biological processes terms (Supplemental Digital Content (Table S1)). Among them, dysfunctional genes were implicated in macrophage activation involved in immune response, T-helper 17 type immune response, and other related biological functions (Fig. 3A). On the other hand, enrichment analysis of KEGG pathway revealed 157 terms (Supplemental Digital Content (Table S1)), including p53 signaling pathway, chronic myeloid leukemia, etc. (Fig. 3B).These signaling pathways played a vital role in the process of vitiligo. Based on GSVA, dysregulation of these signaling pathways was closely related to vitiligo progression (Fig. 3C). Moreover, the signal pathways dysregulated in vitiligo were also explored through GSEA. The results showed that 6 signal pathways were activated incessantly from Non-Lesional to Peri-Lesional and then to Lesional (Fig. 3D–F). In the dataset GSE65127, 4 of the dysregulated signal pathways were verified by an independent dataset GSE75819 (Fig. 3G, H). These results indicated that the signal pathways such as VEGF and mTOR were continuously dysregulated in the pathogenesis of vitiligo.

Figure 3

The enrichment analysis of signal pathways of dysregulated gene in vitiligo. A, The enrichment of biological function of modular dysregulated genes. B, The enrichment of signal pathways of modular genes in KEGG. From blue to red, the significance of the nodes gradually increases. C, Thermogram shows the activation state of signal pathways of KEGG in different groups after processing by GSVA. The red node represents upregulation and the blue node represents downregulation. The common signaling pathways in GSEA for Non-Lesional (D), Peri-Lesional (E), and Lesional (F). G, The signaling pathways in GSEA for Lesional in dataset GSE65127. H, The identical signaling pathways of Lesional in dataset GSE65127 and GSE75819. GSEA = gene set enrichment analysis, GSVA = gene set variation analysis, KEGG = Kyoto Encyclopedia of Genes and Genomes.

Persistent maladjusted genes in the development of vitiligo

In case there was a persistent dysregulation in the signaling pathway, there may be genes with persistent disorders. Consequently, 1197 genes were identified in the 3 groups of DEGs (Fig. 4A). Among them, 269 genes were upregulated constantly from Non-Lesional to Peri-Lesional and then to Lesional, while 82 genes were downregulated (Supplemental Digital Content (Table S2)). Hence, we held that these persistent dysregulated genes may be involved in the occurrence and development of vitiligo. Notably, the expression pattern of persistent maladjusted genes in the modules was displayed by thermogram (Fig. 4B). It was obvious that most of the persistent maladjusted genes were clustered in the module blue. Remarkably, HSPB6 was also a gene that continued to be downregulated in this study. In addition, the continuous dysregulated genes in the module had a tendency of up- or down-regulation from Non-Lesional to Lesional. Significantly, the persistent maladjusted genes were found to participate in T-cell receptor signaling pathway, etc. These signaling pathways functioned as a crucial part in the development of vitiligo. The results revealed that there was a constant alteration of gene expression in the course of vitiligo, which may be related to the severity of vitiligo.

Figure 4

Persistent maladjusted genes in 3 groups of differentially expressed genes. A, Venn diagram of the 3 groups of differentially expressed genes, DEG1 represents Non-Lesional, DEG2 represents Peri-Lesional, and DEG3 represents Lesional. B, Thermogram displays the expression of persistent maladjusted genes in the modules. Red node represents upregulated gene, blue node represents downregulated gene.

Regulation network of persistent maladjusted genes in vitiligo

Transcription regulation indicates that the level of gene expression alters with the alteration of the transcription rate, playing an essential role in transmission of genetic information accurately and diversely.[ To this end, transcription regulation of module genes was explored, which showed that 42 transcription factors had significant regulatory effect on the module genes (Supplemental Digital Content (Table S3)). By screening the regulators of persistent maladjusted genes, it was found that specificity protein 1 (Sp1) was involved in the development of vitiligo by targeting ITGAV to regulate the thyroid hormone signaling pathway. Alternatively, based on drug prediction, 228 drugs were spotted which may have therapeutic effects on module genes (Supplemental Digital Content (Table S4)). Subsequently, focusing on persistent dysregulated genes, transcription factors and drugs were extracted and a map of regulatory network was delineated (Fig. 5A). Furthermore, as expected, the expression of these key genes was verified in GSE75819 dataset (Fig. 5B). The outcomes deciphered that the persistent dysregulated genes affecting the progression of vitiligo were modulated by a number of transcription factors and drugs. Besides, the expressions of hub genes (CACTIN, DCTN1, GPR143, MRPL47, NKTR) in vitiligo patients were verified by qRT-PCR, which were consistent with our findings (Fig. 5C).

Figure 5

Regulation and verification of key persistent maladjusted genes in vitiligo modules. A, The transcription factors and drugs that regulate the persistent maladjusted genes in the modules, as well as the signal pathway network participated by the target genes. Arrows represent transcription factors, diamonds represent drugs, circles represent target genes, triangles represent modules, and squares represent gene-enriched signaling pathways. B, Data in GSE75819 confirmed the differential expression of key persistent maladjusted genes in vitiligo. The red box represents downregulation while the blue box represents upregulation. C, The qRT-PCR results of hub genes in blood samples of vitiligo patients compared with healthy controls,∗ Indicates P < .05; ∗∗∗ Indicates P < .001.

Discussion

The module analysis methods and other related data were applied to study the potential molecular linkages in different process stages of occurrence and development of vitiligo, that was from Non-Lesional, to Peri-Lesional then to Lesional. Under these circumstances, the data from Lesional tissues, Peri-Lesional tissues, and Non-Lesional tissues compared with normal tissues respectively could represent, to a certain extent, the different process stages of occurrence and development of vitligo.

Accordingly, based on WGCNA, the correlation between module and phenotype during the course of vitiligo was identified, including module brown, module blue, and module green, respectively. The module brown was mainly enriched in response to oxygen levels and Wnt pathway. Studies have shown that oxidative stress does play a vital role in causing subsequent autoimmune responses associated with vitiligo.[ Lately, Tang et al[ demonstrated Vitamin D protected human melanocytes against oxidative damage by activation of Wnt/β-catenin signaling. The module blue was mostly enriched in autophagy and immune response. The damage-associated molecular patterns or autoantigens were generated by stressed melanocytes, which initiated innate immunity and adaptive immunity, activating an inflammatory cascade leading to the dysfunction and death of melanocytes.[ The vulnerability of melanocytes to oxidative stress may be induced by dysfunctional autophagy.[ Furthermore, the module green was largely enriched in divalent inorganic cation transport and calcium ion transport. Calreticulin was an omnipresent endoplasmic reticulum protein that regulated translocation of intracellular calcium ions from endoplasmic reticulum lumen to the melanocytes surface under oxidative stress.[ Besides, calreticulin also prompted pro-inflammatory cytokines relating with immune responses, including interleukin-6 and tumor necrosis factor-α.[

Remarkably, for the current study, persistent maladjusted genes from Non-Lesional to Lesional were detected, which may indicate the occurrence and development of vitiligo. Enrichment analysis showed that these genes were significantly enriched in T-cell receptor signaling pathway. Raam et al described in detail that vitiligo was a T-cell-mediated secondary autoimmune disease after treatment of alemtuzumab.[

Combined with gene difference analysis and network topology, potential hub genes were identified. Up to present, the relationships between the 8 hub genes identified in this study and vitiligo have not been reported in the literature. They may be potential biomarkers and therapeutic targets for vitiligo in the future. The hub gene of module brown was NUF2 (NUF2 component of NDC80 kinetochore complex), also known as: cell division cycle associated 1. In line with previous findings, Thang et al[ discovered that knockdown of cell division cycle associated 1 in oral cavity carcinoma cells led to impaired growth and apoptosis in cells. Unfortunately, in the blood samples of vitiligo patients, the expression of NUF2 detected by qRT-PCR was downregulated compared with normal blood samples, which was inconsistent with our findings. However, the expression of NUF2 was verified in an independent dataset GSE75819, in which the sample type was skin tissue. The hub gene of module blue was MRPL47 (mitochondrial ribosomal protein L47). In eukaryotic cells, mitochondria have their own ribosomes (mitoribosomes), which synthesize a small number of proteins, all indispensable for the biogenesis of the oxidative phosphorylation system.[ The hub gene of module green was GPR143 (G protein-coupled receptor 143), which was targeted to melanosomes in pigment cells. It was reported that mutations of GPR143 were connected with X-linked ocular albinism type 1.[ Although skin pigmentation was slightly reduced, subcellular abnormalities such as occurrence of giant melanosomes (macromelanosomes), decrease in melanosome number, and changes of melanosome motility were still observed in epidermal melanocytes.[ The expressions of MRPL47 and GPR143 were confirmed both in the dataset GSE75819 and in the blood samples of vitiligo patients detected by qRT-PCR, which were consistent with our analysis.

Importantly, HSPB6 (heat shock protein family B (small) member 6), as the hub gene of module yellow, was also a persistent maladjusted gene. With mounting evidence accumulated, in response to varied cellular stress or damage, upregulated HSPB6 defended cells against otherwise lethal conditions, such as eliciting smooth muscle relaxation, restraining platelet aggregation, and playing a key regulatory role in cardioprotection and neuroprotection.[ Although the expression level of HSPB6 proteins in normal skin tissues was low, it played a significant role in the detection data in the current study, suggesting that HSPB6 was also an important gene in the pathogenesis of vitiligo. Besides, the ROC curve displayed that HSPB6 as well as other hub genes had potential diagnostic value in dataset GSE65127. However, compared with normal blood samples, the expression of HSPB6 in the blood samples of vitiligo patients was decreased, but there was no statistical difference. By expanding the sample size, it may provide latent biomarker for future research.

In addition to crucial genes, there were also dysregulated signaling pathways that played a vital role in the course of vitiligo, including p53 signaling pathway. Birlea et al emphasized the active participation involved in defined molecular pathways that modulated the balance of melanocytes and keratinocytes between stemness and differentiation states, including p53 and its downstream effectors regulating melanogenesis and Wnt/β-catenin participating in proliferative, migratory, and differentiation in different pigmentation systems.[ Lee et al[ provided that Interferon-inducible T-cell alpha chemoattractant (ITAC) induced the melanocytic migration and hypopigmentation through destabilizing p53 via histone deacetylase 5. Generally, these results were consistent with our integrated analysis.

Moreover, the results presented that transcription factors dominated by SP1 had noteworthy regulatory effects on persistent dysfunctional genes. SP1 is significant to the TATA-less genes, which modulates transcription of manifold target genes related to cell growth, differentiation, apoptosis, and immune response.[ Previous researches unveiled that SP1 controlled the maximal constitutive activity and basal expression of Ah receptor (AHR) gene by binding to the GC-rich motifs.[AHR-129 C > T variant could lead to allele-specific binding of SP1 to AHR promoter and further modify its transcription and downstream effectors in vitiligo.[ Based on persistent maladjusted genes, drug prediction and drug target information suggested that Celecoxib had important regulatory effects on dysfunction modules. Celecoxib is the first cyclo-oxygenase 2-selective inhibitor, which has been approved for treating inflammation and pain for a long time.[ Several drugs required further investigation for the therapy or side effects.

In the present study, WGCNA was performed to analyze the DEGs in vitiligo for the first time. In addition, the data acquired from Lesional skins, Peri-Lesional skins, and Non-Lesional skins compared respectively with normal skins may represent the different process stages of occurrence and development of vitligo. However, the main limitation of the present study was the small sample size and there was little literature to support our findings. Regrettably, at this moment, due to the unavailability of tissue samples, it was difficult for us to detect the protein expression of each hub gene. Therefore, we collected blood samples of vitiligo patients for verification of hub genes. However, heterogeneity between different sample types and individuals may lead to inconsistent results, and future larger cohort validation is needed. Taken together, the hub genes and dysregulated signaling pathways as well as transcription factors and drugs may be verified as potential molecular targets and their roles in vitiligo need to elucidate by further studies.

Conclusions

In general, based on WGCNA and a series of comprehensive analyses, functional modules were identified and persistent maladjusted genes in the course of vitiligo were also detected. Hub genes of each module were uncovered, which may develop to be diagnostic biomarker or therapeutic target of vitiligo. At least, the present study might provide an integrated and in-depth insight for exploring the underlying mechanism of vitiligo occurrence and development and predicting potential treatment methods and mechanism as well.

Author contributions

Conceptualization: Xiaojing Kang, Zixian Lei, Shirong Yu, Yuan Ding, Junqin Liang.

Methodology: Zixian Lei, Shirong Yu, Yilinuer Halifu.

Software: Zixian Lei, Shirong Yu.

Validation: Xiaojing Kang, Zixian Lei, Shirong Yu, Yuan Ding, Junqin Liang.

Formal analysis: Zixian Lei, Shirong Yu, Yuan Ding, Junqin Liang, Yilinuer Halifu, Fang Xiang, Dezhi Zhang.

Investigation: Zixian Lei, Shirong Yu, Yuan Ding, Yilinuer Halifu, Fang Xiang, Dezhi Zhang.

Resources: Zixian Lei, Hongjuan Wang, Wen Hu, Tingting Li, Yunying Wang, Xuelian Zou, Kunjie Zhang.

Data Curation: Zixian Lei, Yilinuer Halifu, Fang Xiang, Dezhi Zhang, Hongjuan Wang, Wen Hu, Tingting Li.

Writing - Original Draft: Zixian Lei, Shirong Yu.

Writing - Review & Editing: Xiaojing Kang, Zixian Lei.

Visualization: Zixian Lei, Yunying Wang, Xuelian Zou, Kunjie Zhang.

Supervision: Xiaojing Kang, Yuan Ding, Junqin Liang.

Project administration: Xiaojing Kang.

Funding acquisition: Xiaojing Kang.

All authors agreed to be accountable for all aspects of this work.