RNA-seq Reveals Dysregulation of Novel Melanocyte Genes upon Oxidative Stress: Implications in Vitiligo Pathogenesis

Abstract

Oxidative stress is known to induce melanocyte death, but the underlying mechanisms are incompletely understood. To identify oxidative stress-induced global gene expression changes in melanocytes, we treated PIG1 melanocytes with H₂O₂ in a dose- and time-dependent manner and performed RNA-seq. This approach allowed us to capture the events occurring early as well as late phase after treatment with H₂O₂. Our bioinformatics analysis identified differentially expressed genes involved in various biological processes of melanocytes which are known to contribute to the vitiligo development, such as apoptosis, autophagy, cell cycle regulation, cell adhesion, immune and inflammatory responses, melanocyte pluripotency, and developmental signaling such as WNT and NOTCH pathways. We uncovered several novel genes that are not previously described to be involved in melanocytic response to stress nor in vitiligo pathogenesis. Quantitative PCR and western blot analysis of selected proteins, performed on PIG1 and primary human epidermal melanocytes, confirmed the RNA-seq data. Interestingly, we discovered an aberrant regulation of several transcription factors that are involved in diabetes, neurological, and psychiatric diseases, all of which are comorbid conditions in patients with vitiligo. Our results may lead to a better understanding of the molecular mechanisms underlying vitiligo pathogenesis and help developing new drug targets for effective treatment.

This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited. The publication of this article was funded by Qatar National Library.

1. Introduction

Reactive oxygen species (ROS) such as hydrogen peroxide (H₂O₂), superoxide anion radical, and hydroxyl radical are generated in the cells endogenously as well as through exposure to extrinsic factors. Under physiological conditions, cells can maintain the intracellular redox homeostasis by scavenging the ROS. However, excessive ROS production disrupts the redox homeostasis and damages the organelles and biomolecules, resulting in the manifestation of a variety of diseases including skin conditions. Thus, ROS can affect diverse biological processes through multiple mechanisms [1].

Skin interfaces with the environment and thus is a major source of ROS. Additionally, ROS are continuously generated during the melanogenesis process in epidermal melanocytes, and this excessive ROS can lead to melanocyte cell death, resulting in skin conditions such as vitiligo [2]. Vitiligo is a progressive skin condition in which functional melanocytes in the epidermis are stressed and selectively destroyed leading to the absence of melanin and a consequent skin depigmentation. Numerous hypotheses about the etiology of vitiligo have been proposed, but it remains unclear what causes damage or death of melanocytes.

There is a compelling evidence that increased production of ROS and their accumulation is one of the major reasons for the death of melanocytes in vitiligo [2]. Very high levels of H₂O₂ have been reported in the epidermis and serum of vitiligo patients [3, 4]. Compared to melanocytes from healthy individuals, those from vitiligo patients showed increased sensitivity and cell death to oxidative stress caused by UVB and chemicals [5]. These observations suggest that melanocytes in vitiligo skin are inherently sensitive to oxidative stress. Thus, identifying and targeting signaling pathways activated by ROS that are responsible for melanocyte death will help prevent the cell damage and probably restore melanocytes in the vitiliginous skin.

Although some of the pathways that underlie H₂O₂-induced transcriptional changes in melanocytes are known, a comprehensive view of stress-induced gene regulation is still elusive. While much of the understanding of ROS-induced signaling in melanocytes is based on single gene/pathway approaches [6, 7], studies on whole-transcriptomic approaches in melanocytes were performed using a single concentration/time point of oxidative stress-inducing agent [8]. Thus, the goal of the present study was to unravel, using RNA-seq, diverse differentially expressed genes (DEGs) in melanocytes that are modulated by H₂O₂-induced oxidative stress in a time- and dose-dependent manner.

2. Results and Discussion

2.1. Generation of Oxidative Stress in Melanocytes Using H₂O₂ and Its Detection

Various batches of isolated primary human epidermal melanocytes (HEM) respond differently to same treatments with various drugs including H₂O₂. These limitations can be avoided by using immortalized human epidermal melanocyte (PIG1) cell line which is a well characterized and most widely used for the fundamental understanding of melanocyte biology [7, 9, 10]. PIG1 cells and primary HEM exhibited dendriticity and extensive DOPA oxidase (tyrosinase) activity and expressed melanocyte-specific marker TYRP-1 and a master regulator of melanocyte development MITF transcription factor, all of which are characteristic features of melanocytes (Figures 1(a) and 1(b)). As compared to untreated controls, H₂O₂ triggered the ROS production, in a concentration-dependent manner as shown by progressive elevation of orange signal in treated cells (Figure 1(c)). Within 3 h after exposure to as low as 100 μM H₂O₂, cells start losing their dendriticity, a first hallmark of cell death process (Figure 1(d)).

[figures omitted; refer to PDF]

2.2. H₂O₂-Induced Cell Death in PIG1 Melanocytes

To choose the right concentration of H₂O₂ and dosing time points for RNA-seq experiments, we measured the cytotoxicity of H₂O₂ in melanocytes by exposing them to increasing concentrations of H₂O₂ for various exposure times. Cell viability was found to decrease in an inversely correlated manner to both dose and time of H₂O₂ exposure (Figure 1(e)). Across all concentrations of H₂O₂, cell death was minimal at time point 6 h. While a significant decrease in cell viability to 25% (24 h) and 40% (48 h) was observed in the presence of 250 μM H₂O₂, it was reduced by 50% and 90% with 350 μM H₂O₂ at these respective time points. Beyond these concentrations and time points, no cell viability was observed. Thus, to determine changes in gene expression using RNA-seq, wherein all stages of death process (early, mid, and late stages) could be represented, we have chosen 100 μM (nonlethal dose, generally seen in healthy cells), 250 μM (moderately lethal dose), and 350 μM (highly lethal) of H₂O₂ and time points of 6 h, 24 h, and 48 h and performed with biological triplicates (Figure 1(f)).

2.3. Gene Expression Pattern

MDS plot shows a clear clustering of controls and H₂O₂-treated samples (Figure 2(a)). The volcano plot (Figure 2(b)) shows an overall representation of RNA-seq results obtained with cells treated with 250 μM H₂O₂ for 48 h. Similar pattern is observed with all other treatments (not shown).

[figures omitted; refer to PDF]

2.4. H₂O₂ Alters Gene Expression in a Time- and Dose-Dependent Manner in Melanocytes

H₂O₂ altered the expression changes in a time- and dose-dependent manner. We noticed a general trend that the number of genes induced by H₂O₂ exceeded the number of suppressed ones across time and concentrations (Figure 2(c)) (hypergeometric $p$ value < 0.05 for DEGs between time points and concentrations). We noticed that there are significant overlapping DEGs between edgeR, DESeq2, and limma methods (Sup Fig-1). Figure 3(a) shows the hierarchical clustering of top 100 DEGs between controls and 250 μM H₂O₂ treatment, and Sup. Fig-2 shows hierarchical clustering at other concentrations. The list of top 25 upregulated and downregulated genes (250 μM, 48 h) is shown in Tables 1 and 2, respectively. While top upregulated genes are involved in calcium signaling, GPCR signaling, and growth and differentiation, top downregulated genes are involved in cytoskeleton, adhesion, immune mediators, and oxidation-reduction pathways.

[figures omitted; refer to PDF]

Table 1

List of top 25 most upregulated genes by treatment with 250 μM H₂O₂ for 48 h. These genes passed the $\text{FDR}\text{threshold}<0.05$ and abs.log₂FC of 2.

Gene_ID	Gene symbol	Gene name	Log2 fold change
ENSG00000120903	CHRNA2	Cholinergic receptor nicotinic alpha 2 subunit	7.1
ENSG00000100095	SEZ6L	Seizure related 6 homolog like	6.5
ENSG00000198576	ARC	Activity regulated cytoskeleton associated protein	6.2
ENSG00000124253	PCK1	Phosphoenolpyruvate carboxykinase 1	6.2
ENSG00000137270	GCM1	Glial cells missing homolog 1	6.1
ENSG00000184647	PRSS55	Serine protease 55	5.9
ENSG00000226364	AC102948.2	No description	5.9
ENSG00000182870	GALNT9	Polypeptide N-acetylgalactosaminyltransferase 9	5.7
ENSG00000087237	CETP	Cholesteryl ester transfer protein	5.3
ENSG00000142233	NTN5	Netrin 5	5.3
ENSG00000168081	PNOC	Prepronociceptin	5.0
ENSG00000172421	EFCAB3	EF-hand calcium binding domain 3	4.8
ENSG00000172137	CALB2	Calbindin 2	4.8
ENSG00000182896	TMEM95	Transmembrane protein 95	4.7
ENSG00000178623	GPR35	G protein-coupled receptor 35	4.7
ENSG00000185739	SRL	Sarcalumenin	4.7
ENSG00000178297	TMPRSS9	Transmembrane serine protease 9	4.7
ENSG00000136099	PCDH8	Protocadherin 8	4.6
ENSG00000225899	FRG2B	FSHD region gene 2 family member B	4.6
ENSG00000167772	ANGPTL4	Angiopoietin like 4	4.5
ENSG00000168746	C20orf62	Uncharacterized protein C20orf62	4.5
ENSG00000164266	SPINK1	Serine peptidase inhibitor, Kazal type 1	4.4
ENSG00000188886	ASTL	Astacin-like metalloendopeptidase	4.3
ENSG00000156466	GDF6	Growth differentiation factor 6	4.3
ENSG00000115353	TACR1	Tachykinin receptor 1	4.3

Table 2

List of top 25 most downregulated genes by treatment with 250 μM H₂O₂ for 48 h. These genes passed the $\text{FDR}\text{threshold}<0.05$ and abs.log₂FC of 2.

Gene ID	Gene symbol	Gene name	Log2 fold change
ENSG00000205420	KRT6A	Keratin 6A	-8.6
ENSG00000171345	KRT19	Keratin 19	-8.3
ENSG00000163220	S100A9	S100 calcium-binding protein A9	-7.9
ENSG00000147256	ARHGAP36	Rho GTPase activating protein 36	-7.6
ENSG00000148346	LCN2	Lipocalin 2	-7.6
ENSG00000086548	CEACAM6	Carcinoembryonic antigen related cellAdhesion molecule 6	-7.0
ENSG00000168542	COL3A1	Collagen type III alpha 1 chain	-7.0
ENSG00000057149	SERPINB3	Serpin family B member 3	-6.7
ENSG00000087916	SLC6A14	Solute carrier family 6 member 14	-6.6
ENSG00000124107	SLPI	Secretory leukocyte peptidase inhibitor	-6.6
ENSG00000173432	SAA1	Serum amyloid A1	-6.4
ENSG00000186081	KRT5	Keratin 5	-6.4
ENSG00000124102	PI3	Protease inhibitor 3, skin-derived	-6.3
ENSG00000164107	HAND2	Heart and neural crest derivatives expressed 2	-6.3
ENSG00000181541	MAB21L2	Mab-21 like 2	-6.1
ENSG00000125810	CD93	CD93 molecule	-6.0
ENSG00000127954	STEAP4	STEAP4 metalloreductase	-6.0
ENSG00000158578	ALAS2	5${}^{\prime }$-Aminolevulinate synthase 2	-5.9
ENSG00000162896	PIGR	Polymeric immunoglobulin receptor	-5.9
ENSG00000016082	ISL1	ISL1 transcription factor, LIM/homeodomain	-5.9
ENSG00000186832	KRT16	Keratin 16	-5.9
ENSG00000125999	BPIFB1	BPI fold containing family B member 1	-5.9
ENSG00000129538	RNASE1	Ribonuclease 1	-5.8
ENSG00000171401	KRT13	Keratin 13	-5.8
ENSG00000166482	MFAP4	Microfibril-associated protein 4	-5.8

2.5. GO Enrichment and Pathway Analysis

To better understand the functions of DEGs, we performed GO enrichment analysis using DAVID. At the level of molecular functions, the top 3 enrichment items included the ion binding, endopeptidase activity, and transcriptional activator activity (Figure 3(b) and Sup. Fig-3), whereas cell adhesion, immune response, and extracellular matrix organization were the most enriched in biological process. At the level of cellular component, the integral component of membranes and plasma membrane was top enriched. Pathways related to neuronal signaling, immune functions and inflammatory response, apoptosis/cell survival such as PI3K and MAPK pathways, cell adhesion, and migration are among the most relevant enriched ones (Figure 3(c), Sup. Fig-4). Similar results have been reported in previous studies that compared the differences in expression pattern in healthy vs. vitiligo skin [11].

Analysis of DEGs shows that the expression of certain genes known to be involved in signaling mechanisms contributing to various melanocytic processes is quite significant across time points in H₂O₂-treated cells. These were classified by their potential roles in cell death, cell cycle, melanogenesis, WNT signaling, etc. (Table 3). Since some of these genes are known to involve in more than one biological process, these were classified into several groups. Furthermore, the most interesting part of our study, we found several genes that are not previously described to be involved in melanocyte death or in vitiligo pathogenesis (Table 3). A detailed account of these novel genes and their possible role in melanocyte survival has been discussed in following sections.

Table 3

Functional classes of selected differentially expressed genes with 250 μM H₂O₂ for 48 h. These genes passed the $\text{FDR}\text{threshold}<0.05$ and abs.log₂FC of 2.

Putative function	Genes	Log2 fold change
Apoptosis	GADD45B	+2.6
	MACC1	-5.2
	HGF	-2.6
	GATA4	-2.0
	TP53AIP1	+1.9
	FOS	+2.6
	JUN	+2.1
	JUNB	+1.2
	FOSB	+2.6
	CARD9	+1.9
	ANXA3	-2.7
	IRF8	+3.7
	TCF4	-5.5
	HSPA1A	+2.0
	HSPA1B	+2.2
	CLU	+3.2
	TRAF1	+3.0
	TNFRSF4	+1.3
	TNFRSF12A	+1.5
	TNFRSF13C	+1.8
	GZMM	+3.5
Autophagy	GATA4	-2.0
	ATG9B	+1.7
Melanogenesis	TYRP1	-1.6
	PMEL	-1.4
	PLP1	-1.8
	HGF	-3.8
Cell cycle	PPP2R2C	+1.9
	CCND2	-3.0
	CDKN1A/p21	+0.7
WNT signaling	DKK1	+1.4
	TCF4	-5.5
	CTNNB1 (b-catenin)	-0.6
Adhesion	MACC1	-5.2
	CDH1/E cadherin	-1.8
	PTGS2	-1.4
Stress response	HSPA1A	+2.0
	HSPA1B	+2.2
Immunity	IL1B	-4.7
	IL6	+2.1
	IL17F	+3.3
	CXCL10	-5.0
	CCL2	-4.0
	CLU	+3.2
	TNFRSF4	+1.3
	TNFRSF12A	+1.5
	TNFRSF13C	+1.8
	TNFRSF25	+1.3
	IRF8	+3.7
	OASL	-1.8

- = downregulated; + = upregulated genes.

2.6. Validation of RNA-seq Data

We performed quantitative RT-PCR of selected genes on PIG1 melanocytes as well as primary HEM treated with 250 μM H₂O₂ for 48 h, using GAPDH as a reference gene. At these conditions, the cell death in PIG1 and primary HEM cells was 36% and 31%, respectively (Figure 4(a)). We found an excellent agreement between the RNA-seq data and RT-PCR (Figure 4(b)). Similar results were observed when β-actin was used as a reference gene (Sup. Fig-5). Furthermore, we checked the expression of selected genes, CRAD9, TCF4, FOS, KIT, CLU, and MACC1, at protein level using the cell lysates of control and H₂O₂-treated PIG1 cells and primary HEM by western blotting. The expression patterns of these proteins (Figure 4(c)) confirm that there exists a good concordance between RNA-seq, RT-PCR, and western analysis, thus enhancing our confidence that melanocyte death is mediated by the differential molecular abnormalities of these proteins.

[figures omitted; refer to PDF]

2.7. Cell Death Mechanisms

2.7.1. Apoptosis

Apoptosis has been reported to be the main mechanism of melanocyte destruction by ROS in vitiligo [12–16]. Furthermore, previous GWAS studies identified that several candidate loci associated with vitiligo pathogenesis are regulators of apoptosis [17]. Apoptosis can occur through two distinct molecular pathways: intrinsic or extrinsic pathways.

(1) Intrinsic Pathway. Lower expression of antiapoptotic proteins and elevated expression of proapoptotic proteins such as BAX and p53 have been noticed in vitiliginous skin as compared to normal skin [18]. For the first time, our RNA-seq results suggest that the regulation of apoptosis in stress-induced melanocytes is more complex than previously reported. Thus, while several proapoptotic proteins such as BAX, BAD, BIM, BID, BIK, BOK, HRK, NOXA, and PUMA were found to be upregulated in response to H₂O₂-induced stress, we found that the expression of antiapoptotic members such as BCL11A, BCL11B, A1, and API5 was suppressed (Figures 4(d) and 5). While differential expression of each of these members has a modest effect, simultaneous elevation of multiple proapoptotic genes and downregulation of several antiapoptotic genes will tilt the balance towards apoptosis in response to stress.

[figure omitted; refer to PDF]

(2) Extrinsic Apoptotic Pathway. As far as extrinsic pathway is concerned, the members of the tumor necrosis factor receptor superfamily (TNFRSF) bind to death ligands TNFs. They are primarily involved in diverse biological processes such as immune homeostasis, execution of immune responses, inflammation, stimulation of apoptosis, and proliferation [19]. The most interesting observation from our study is that several members of TNFRSF such as TNFRSF-1B, 4, 8, 9, 10A, 11B, 12A, 13C, and 25 are significantly upregulated to various extent after treatment with H₂O₂ (Table 3, Sup. Table-1). The TNF-α, the ligand that binds to TNFRSF, has been indeed shown to accumulate in the skin and serum of vitiligo patients [20]. The overexpression of TNFRSF members may have a huge effect on cells. While one way, they can help execute immune responses against oxidative stress, on the other way, they activate melanocyte cell death.

2.7.2. Autophagy

In addition to apoptosis, H₂O₂ also induced the expression of genes involved in autophagy. Of these, downregulation of a zinc finger TF, GATA4, is worth a mention. While silencing of GATA4 can trigger autophagy and apoptosis, overexpression of GATA4 elevated the gene expression of the survival proteins and suppressed the expression of other autophagy-related genes [21]. Suppression of GATA4 by H₂O₂ as seen in our study, an observation consistent with a previous report showing the downregulation of GATA3 in vitiligo melanocytes [22], may likely be responsible, at least partially, for the autophagic effects of H₂O₂. ATG9B, whose up expression was more prominent at 24 h after treatment with H₂O₂, is also well known to participate in autophagy [23]. Other autophagy genes either unchanged or downregulated may suggest that autophagy in H₂O₂-stressed melanocytes preferentially depend on GATA4 and ATG9B.

2.7.3. Melanogenesis

In addition to the melanocyte death, abnormal melanogenesis is thought to contribute to the vitiligo pathogenesis. Consistent with previous studies [24, 25], we noticed a downregulation of several genes involved in pigmentation process, such as TYRP1, PMEL, MLANA, DCT, and PLP1, whose underexpression was more prominent at 48 h after treatment with H₂O₂ (Table 3), suggesting that these genes are aberrantly regulated by oxidative stress and play a role in disease pathology.

2.7.4. Other Novel Cell Death Signaling Pathways

Besides classical BCL2 family members, numerous other proteins are known to control cell survival or death either directly or indirectly. For the first time, our study identified several such proteins and therefore it is worth discussing the most significant of them and their possible implications in melanocyte biology.

It has been shown that metastasis associated in colon cancer protein 1 (MACC1) can promote cell growth when overexpressed and promoted apoptosis in both in vitro and in vivo when underexpressed [26]. MACC1 activates the HGF/c-MET pathway, culminating in aberrant activation of multiple cellular responses such as proliferation, cell morphogenesis, migration, and breakdown of the extracellular matrix by altering Ras/MAPK and PI3K/Akt signaling pathways [27]. A key finding of the present study was that H₂O₂-induced oxidative stress significantly reduced the expression of both MACC1 and HGF in melanocytes (Table 3). This probably leads to the repression of Ras/MAPK and PI3K/Akt survival pathways, resulting in suppression of cell proliferation and induction of apoptosis (Figure 5). Further experiments are in progress to identify the underlying mechanism of MACC1/HGF-mediated apoptosis in melanocytes.

Growth Arrest and DNA Damage Inducible (GADD) family proteins are implicated in cell cycle arrest, apoptosis, innate immunity, and maintenance of genomic stability [28]. The transcription of GADD family proteins is induced by apoptotic cytokines and genotoxic stress. These proteins mediate apoptosis by activating the p38/JNK pathway. We discovered a concomitant upregulation of both GADD45B and p38 MAPK in melanocytes treated with H₂O₂ (Table 3), suggesting that at least part of the effects of GADD45 proteins on cell growth and apoptosis are mediated by activation of the p38 pathway.

Caspase Activation and Recruitment Domain (CARD) are protein interaction motifs found in a variety of proteins such as CARD9, CARD11, and CARD14. They are known to participate in activation or suppression of CARD containing members of the caspase family and play a pivotal regulatory role in cell apoptosis and inflammation [29]. We found that all these proteins are upregulated more significantly at 24 h after treatment with 250 μM H₂O₂ (Sup. Table-1). This is consistent with a previous observation showing that CARD11 is upregulated in the Smyth line of chicken, which is an excellent avian model for human autoimmune vitiligo [30]. Similarly, variations in CARD7/NALP1 gene are associated with the development of generalized vitiligo, and its overexpression was demonstrated to induce apoptosis [31]. While NALP1 expression is only modestly increased, the expression of NALP3, another member of this family with similar function, is substantially increased in melanocytes exposed to H₂O₂ at 24 h (Sup. Table-1).

The expression of Tumor Protein p53-Regulated Apoptosis-Inducing Protein 1 (TP53AIP1) gene is inducible by p53 and is thought to play an important role in mediating p53-dependent apoptosis. An increased level of proapoptotic protein p53 was found in the lesional skin compared to perilesional or nonlesional areas in vitiligo patients [12]. In our study, a higher expression of TP53AIP1 was apparent at 48 h after treatment with H₂O₂, pointing to the fact that in addition to the direct involvement of aberrantly expressed BCL2 family proteins, the apoptosis in melanocytes could be induced by p53-dependent pathways.

2.7.5. cKIT Expression, Heat Shock Response, and FOS/JUN Signaling

A modest decrease in the expression of c-KIT receptors in melanocytes treated with H₂O₂ was observed in our study. This may contribute to the melanocyte growth arrest and death, consistent with previous reports showing the downregulation of cKIT in the skin of vitiligo [24]. Furthermore, our results showed that stressed melanocytes expressed higher levels of HSPA1A/HSP70-1 and HSPA1B/HSP70-2 and FOS/JUN/p38 MAPK proteins, all of which can contribute to cell death (Figure 5).

The proapoptotic TNF signals are blocked by proteins that are induced by NF-κB such as TNFR-associated factor 1 (TRAF1) [32]. Our results suggest that elevation of TRAF1 acts as an antiapoptotic mechanism to prevent death of melanocytes in response to stress. However, the enhanced expression of NFKB Inhibitor Zeta (NFKBIZ), an antagonist of NFKB, along with other proapoptotic proteins present in stressed cells, as observed in our study, possibly negate the protection offered by TRAF and ultimately ensues cell apoptosis.

2.7.6. Cell Cycle Regulating Proteins

D-type cyclins (cyclin D1, D2, and D3) are well known to play critical roles in cell cycle progression by interacting with cyclin-dependent kinases, such as CDK4 and CDK6. High expression of cyclins has been detected in several tumors. Repression of cyclin D2 (CCND2) expression, and concomitant upregulation of PPP2R2C, that is implicated in the negative control of cell growth and division [33], as we observed in cells treated with 250 μM H₂O₂ for 48 h, could result in the G1 arrest and subsequent growth retardation. The relatively unchanged or decreased expression of other cyclins may suggest that the arrest of the cell cycle progression in stressed cells preferentially depend on CCND2. Together, our results suggest that H₂O₂ induces cell cycle arrest by regulating the expressions of cell cycle-related proteins.

2.7.7. Developmental/Pluripotency Pathways

Our RNA-seq analysis also revealed the modulation of proteins involved in developmental/pluripotency signaling, such as WNT, and NOTCH pathways. WNT5A, 5B, 6, 10B, and 11 and DKK1 are among the most differentially regulated ones by treatment with H₂O₂. WNT signaling is critical for the development of melanocyte and was shown to be affected in vitiligo skin [24]. Upregulation of WNT pathway inhibitor, DKK1 in our study was most impressive. This is because while incubation of melanocytes with DKK1 resulted in increased apoptosis and reduced pigmentation, overexpression of DKK1 was previously shown to suppress melanocyte function and growth by reducing the expression levels of MITF, DCT, tyrosinase, and PMEL [34]. This is consistent with our findings as we demonstrated that the expression levels of all these genes were indeed repressed. Thus, in addition to activating pathways that are directly involved in melanocyte death, H₂O₂ is also inhibiting the regeneration of melanocytes by altering the signaling mechanisms involved in melanocyte development.

Transcription factor 4 (TCF4) can control the nuclear response to Wnt/β-catenin signaling and also functioning of immune system cells, neurons, and melanocytes. Knockdown of TCF4 has been shown to induce cell cycle arrest and apoptosis [35]. Since significant suppression of TCF4 by H₂O₂ was evident from our study, it might be possible that TCF4-mediated cell cycle arrest and apoptosis form another layer of regulatory circuitry for the destruction of melanocyte under oxidative stress or in vitiligo condition. Other proteins that we found differentially regulated by H₂O₂ include ZEB1, SNAI1, MYT1L, SOX21, KLF4, and GRHL2, all of which are well known to control epithelial to mesenchymal transition, development, or pluripotency.

2.7.8. Inflammation and Immune Signaling

Interleukin-17 is a family of six cytokines that includes IL-17A through IL-17F, which are well characterized for their roles in immune modulation. Increased expression of IL17 in the serum and skin of vitiligo patients has been reported [36]. In our present work, we found that among IL17 family cytokines, IL17F is substantially upregulated in melanocytes treated with H₂O₂. This upregulation can antagonize melanogenesis and also promote melanocyte death by downregulating BCL2 family proteins, an observation consistent with previous reports [37]. Similarly, serum concentration of IL6 is elevated in vitiligo patients [38], and our results confirmed IL6 overexpression at a cellular level. IL6 can induce the expression of cell adhesion molecules in melanocytes, thereby facilitating the interaction of melanocytes with immune cells and possibly induces B-cell activation, increasing the autoantibody production and subsequent damage of melanocytes [39]. In addition to IL6, our results show that several other proteins involved in immunity and inflammation, such as CXCL17 (-4.7), IL1A (-4.8), IL1B (-4.7), IL1RN (-3.5), and OASL (-1.8), have been found to be aberrantly expressed in stressed melanocytes [40].

2.7.9. Transcription Factors

Transcription factors are critical for the transcriptional regulation of gene expression and play a key role in many biological processes. Therefore, we analyzed our RNA-seq data to obtain aberrantly regulated TFs in response to oxidative stress. We identified 30 upregulated and 18 downregulated TFs (Table 4). Some of these differentially regulated TFs play a role in immune system, pluripotency, differentiation, development, and cell death. A novel finding from our study is the aberrant expression in stressed melanocytes of TFs that play a role in diabetes, neurological, and psychiatric conditions, all of which are potential comorbid diseases in patients with vitiligo [41–43], but no role of these TFs has been previously reported in melanocytes.

Table 4

Differentially regulated transcription factors between control and 250 μM H₂O₂ for 48 h and their putative functions. These genes passed the $\text{FDR}\text{threshold}<0.05$ and abs.log₂FC of 2.

Putative function	Genes	Log2 fold change
Neuronal/brain development/psychiatric/mental retardation	POU3F1	+4.0
	MYT1L	+3.8
	DLX3	+3.7
	PITX3	+3.4
	FOXJ1	+3.3
	NHLH1	+2.9
	EGR3	+2.8
	ZNF804B	+2.2
	SOX21	+2.2
	GBX1	+2.0
	BARHL2	+2.0
	TCF4	-5.5
	ZFHX4	-3.4
	MEIS2	-3.1
	SETBP1	-2.6
	TFAP2B	-2.2
Immune system	IRF8	+3.7
	FOXN1	+3.5
	TBX21	+2.7
	IKZF1	-5.3
	IRF6	-2.1
Pluripotency/EMT	MYT1L	+3.8
	SNAI1	+3.6
	SOX21	+2.2
	KLF4	+2.0
	ZEB1	-5.5
	TCF4	-5.5
	GRHL2	-2.3
Growth/proliferation/differentiation/development	FOXN1	+3.5
	SP6	+3.5
	NHLH1	+2.9
	FOXL1	+2.8
	EGR3	+2.8
	FOSB	+2.6
	EGR1	+2.5
	SOX21	+2.2
	OVOL1	+2.1
	EHF	-5.6
Insulin/diabetes	MAFA	+3.7
	INSM2	+2.2
	BARHL2	+2.0
	ISL1	-5.9
Cardiovascular	GATA5	+3.0
	HAND2	-6.2
	MEIS2	-3.1
Metabolism	FOXL1	+2.8
Skin development	OVOL1	+2.1
	KLF4	+2.0
Cell death/survival	SNAI1	+3.6
	CSRNP3	-2.8
	IRF6	-2.1

- = downregulated; + = upregulated genes.

Our findings shed light for the first time on this comorbid relationship: the TFs differentially expressed in stressed melanocytes might be the link between vitiligo and potentially associated conditions. This can be partially explained by the fact that although melanocytes are mainly found in the skin, they are also recognized in other parts of the body, including the eyes, ears, heart, and central nervous system where they are thought to have different roles from that played in the skin [42], but they may react in a similar way to stress. Another explanation is that the skin and brain share the same embryonic ectodermal origin, influenced by same hormones and neurotransmitters, and thus may have a similar behavior to stress. Additionally, nerve, pancreatic, and cardiac cells, like melanocytes, are vulnerable to ROS, and overproduction of ROS is shown to cause several diseases such as diabetes, rheumatoid arthritis, cardiovascular diseases, stroke, cancer, and other degenerative diseases in humans. It could be postulated that melanocytes, nerve, or pancreatic cells may respond to the oxidative stress by activation of common set of TFs, leading to the regulation of certain genes responsible for the occurrence of comorbid diseases. For example, myelin transcription factor 1-like (MYT1L) is known to be specifically expressed in the brain, and its dysregulation was shown to cause intellectual disability like autism spectrum disorders [44]. Similarly, ISL1 (ISL LIM homeobox 1) plays an important role in regulating insulin gene expression, as well as in the development of pancreatic cell lineages and motor neuron generation [45, 46]. Mutations in ISL1 have been associated with maturity onset diabetes of the young and type-2 diabetes [47]. Our results show that both these TFs are differentially regulated in stressed melanocytes, suggesting that they have a role not only in nerve/pancreatic cells as previously reported but also in melanocytes. A functional cross talk between the melanocyte, nervous, and immune systems can be another explanation for the occurrence of comorbid conditions, but further investigation is needed to clarify a potentially shared etiopathogenesis.

The most striking disadvantage of profiling the vitiligo skin being the lack of melanocytes in lesional skin, thus the actual pathways operated in melanocytes that are responsible for their disappearance remain uncovered. This study tries to fulfill the knowledge gaps encountered in previous studies and to identify novel genes or pathways that play a profound role in vitiligo pathogenesis which may have been missed in earlier studies.

Our study is not without limitations. It is always ideal to study gene expression using 3D skin models, which takes into account the effects of microenvironment. However, various technical challenges such as difficulty in quantifying melanocyte death in a mixed population of cells and difficulty in sorting the stressed melanocytes that are undergoing various cell death processes from a 3D skin model prohibit the use of a 3D model.

In summary, our RNA-seq approach identified potential novel melanocyte genes induced by oxidative stress in a dose-and time-dependent fashion. Novel findings from our study include the notable changes in the expression of diverse genes known to play a role in death and several functions in other cell types, but their role in melanocyte biology or death was not previously described. Thus, basing on the current study, it is reasonable to hypothesize that ROS-induced melanocyte damage is regulated by a complex network of diverse, fail-proof, multilayered signaling mechanisms (Figures 4 and 5). Thus, a multipronged approach is needed to effectively counter ROS effects and to prevent melanocyte loss in conditions like vitiligo. Identification of these novel genes provides an additional clue to vitiligo pathogenesis and opens new avenues for further investigation in melanocyte biology. Further studies are needed to unveil the precise function of these genes, which may help develop new drug targets in conditions associated with melanocyte stress such as in vitiligo.

3. Materials and Methods

3.1. Cell Lines and Reagents

PIG1 melanocyte cell line was a gift from Caroline Le Poole, Northwestern University, Chicago, Illinois, USA. Primary HEM, melanocyte culture medium M254, Human Melanocyte Growth Supplement (HMGS), keratinocyte medium M154, and Human Keratinocyte Growth Supplement (HKGS) were obtained from Thermo Fisher Scientific, USA. Keratinocyte cell line HaCaT was from AddexBio (USA). All antibodies and chemicals, unless specified, were obtained from Cell Signaling Technology (Beverly, MA, USA) and Sigma (Milwaukee, WI, USA), respectively.

3.2. Cell Culture and L-DOPA Staining

PIG1 and primary HEM were cultured in M254 medium supplemented with HMGS, at 37°C in a humidified incubator of 5% CO₂. HaCaT cells were cultured in M154 medium supplemented with HKGS. L-DOPA staining of melanocytes was performed as described previously [48].

3.3. Cell Viability Assay and RNA Isolation

Cells were cultured in 10 cm BioCoat tissue culture plates (Corning) until 70% confluence. Then, medium was changed to supplement-free M254 medium and treated with various doses of H₂O₂ for various time points. At the end of incubation, both dead and live cells were collected, and cell viability was measured by trypan blue dye exclusion assay. Remaining cells were used to isolate DNA and RNA using the Total DNA/RNA isolation kit (Qiagen). All cytotoxicity data was representation of three independent experiments in triplicates. The quality and quantity of RNA were measured by OD A260/A280 by NanoDrop.

3.4. Measurement of ROS

ROS production was assessed as per manufacturer’s instructions. In brief, cells were incubated with 10 μM of CellROX Orange dye (Molecular Probes) for 30 min at 37°C, followed by washing twice with PBS. On contact with ROS, the fluorescein is oxidized and emits an orange fluorescence. Cells were viewed and imaged using EVOS fluorescence microscopy.

3.5. RT-PCR and Western Blotting

About 2 μg of RNA was reverse transcribed to generate cDNA using random primers and Superscript III reverse transcriptase (Invitrogen). Quantitative RT-PCR reactions were carried out in triplicates on a QuantStudio 12K Flex Real Time PCR machine. The relative expression of each gene was calculated using the 2-DDCT method with GAPDH as reference. The primers were obtained from Integrated DNA Technologies, and their sequences can be obtained upon request. Fifty microgram of total protein was loaded on 10% SDS-PAGE gel, and standard western blotting protocol was used to detect proteins.

3.6. RNA-seq Data Analysis

Raw paired-end (PE) reads were adapter-trimmed using Trimmomatic [49]. To maximize the number of reads on mRNA, we filtered raw reads against the Human rRNA databases using the SortMeRNA tool [50]. Then, reads were filtered for vector sequences by mapping to NCBI UniVec database using SeqyClean [51]. Then, the reads were aligned to the GRCh37 reference genome using the STAR aligner [52] with default parameters for PE reads, and gene counts were generated using the Subread package tool featureCounts [53]. The limma/voom R Bioconductor packages were used for normalization of read counts and identification of DEGs between treatment and control groups [54]. GO enrichment and pathway analysis of DEGs were performed using the DAVID bioinformatics tool [55]. To identify DEGs, we set a threshold of absolute log₂ fold $\text{change}\left({\log }_2\text{FC}\right)>2.0$ and $\text{FDR}<0.05$. Furthermore, we used different platforms like edgeR and DEseq2 methods to confirm DEGs. We also statistically compared the DEGs between different contrast groups (time points for each concentration) using hypergeometric tests.

3.7. Statistical Analysis

Unless indicated otherwise, data represent the results for assays performed in triplicate, with error bars to S.D.

Acknowledgments

We thank Dr. Caroline Le Poole, Professor of Dermatology, Northwestern University, Chicago, Illinois, USA, for providing PIG1 melanocytes. This study is funded by Sidra Medicine, Doha, Qatar. The publication of this article was funded by the Qatar National Library.

References

[1] M. Schieber, N. S. Chandel, "ROS function in redox signaling and oxidative stress," Current Biology, vol. 24 no. 10, pp. R453-R462, DOI: 10.1016/j.cub.2014.03.034, 2014.

[2] L. Denat, A. L. Kadekaro, L. Marrot, S. A. Leachman, Z. A. Abdel-Malek, "Melanocytes as instigators and victims of oxidative stress," The Journal of Investigative Dermatology, vol. 134 no. 6, pp. 1512-1518, DOI: 10.1038/jid.2014.65, 2014.

[3] A. Jain, J. Mal, V. Mehndiratta, R. Chander, S. K. Patra, "Study of oxidative stress in vitiligo," Indian Journal of Clinical Biochemistry, vol. 26 no. 1, pp. 78-81, DOI: 10.1007/s12291-010-0045-7, 2011.

[4] K. U. Schallreuter, J. Moore, J. M. Wood, W. D. Beazley, D. C. Gaze, D. J. Tobin, H. S. Marshall, A. Panske, E. Panzig, N. A. Hibberts, "In vivo and in vitro evidence for hydrogen peroxide (H 2 O 2 ) accumulation in the epidermis of patients with vitiligo and its successful removal by a UVB-activated pseudocatalase," The Journal of Investigative Dermatology Symposium Proceedings, vol. 4 no. 1, pp. 91-96, 1999.

[5] K. Jimbow, H. Chen, J. S. Park, P. D. Thomas, "Increased sensitivity of melanocytes to oxidative stress and abnormal expression of tyrosinase-related protein in vitiligo," The British Journal of Dermatology, vol. 144 no. 1, pp. 55-65, DOI: 10.1046/j.1365-2133.2001.03952.x, 2001.

[6] O. A. Arowojolu, S. J. Orlow, N. Elbuluk, P. Manga, "The nuclear factor (erythroid-derived 2)-like 2 (NRF2) antioxidant response promotes melanocyte viability and reduces toxicity of the vitiligo‐inducing phenol monobenzone," Experimental Dermatology, vol. 26 no. 7, pp. 637-644, DOI: 10.1111/exd.13350, 2017.

[7] J. A. Mosenson, K. Flood, J. Klarquist, J. M. Eby, A. Koshoffer, R. E. Boissy, A. Overbeck, R. C. Tung, I. C. le Poole, "Preferential secretion of inducible HSP70 by vitiligo melanocytes under stress," Pigment Cell & Melanoma Research, vol. 27 no. 2, pp. 209-220, DOI: 10.1111/pcmr.12208, 2014.

[8] G. Yang, G. Zhang, M. R. Pittelkow, M. Ramoni, H. Tsao, "Expression profiling of UVB response in melanocytes identifies a set of p53-target genes," The Journal of Investigative Dermatology, vol. 126 no. 11, pp. 2490-2506, DOI: 10.1038/sj.jid.5700470, 2006.

[9] I. C. Le Poole, F. M. van den Berg, R. M. van den Wijngaard, D. A. Galloway, P. J. van Amstel, A. A. Buffing, H. L. Smits, W. Westerhof, P. K. Das, "Generation of a human melanocyte cell line by introduction of HPV16 E6 and E7 genes," In Vitro Cellular & Developmental Biology-Animal, vol. 33 no. 1, article 21, pp. 42-49, DOI: 10.1007/s11626-997-0021-6, 1997.

[10] Z. Zhou, C. Y. Li, K. Li, T. Wang, B. Zhang, T. W. Gao, "Decreased methionine sulphoxide reductase A expression renders melanocytes more sensitive to oxidative stress: a possible cause for melanocyte loss in vitiligo," The British Journal of Dermatology, vol. 161 no. 3, pp. 504-509, DOI: 10.1111/j.1365-2133.2009.09288.x, 2009.

[11] P. Wang, Y. Li, H. Nie, X. Zhang, Q. Shao, X. Hou, W. Xu, W. Hong, A. Xu, "The changes of gene expression profiling between segmental vitiligo, generalized vitiligo and healthy individual," Journal of Dermatological Science, vol. 84 no. 1, pp. 40-49, DOI: 10.1016/j.jdermsci.2016.07.006, 2016.

[12] A. M. Abdel-Aal, M. A. Kasem, A. H. Abdel-Rahman, "Evaluation of the role of apoptosis in vitiligo: immunohistochemical expression of P53, Bcl-2 and MART-1 antigens," The Egyptian Journal of Hospital Medicine, vol. 8, 2002.

[13] C. L. Huang, J. J. Nordlund, R. Boissy, "Vitiligo: a manifestation of apoptosis?," American Journal of Clinical Dermatology, vol. 3 no. 5, pp. 301-308, DOI: 10.2165/00128071-200203050-00001, 2002.

[14] R. Kumar, D. Parsad, "Melanocytorrhagy and apoptosis in vitiligo: connecting jigsaw pieces," Indian Journal of Dermatology, Venereology and Leprology, vol. 78 no. 1, pp. 19-23, DOI: 10.4103/0378-6323.90942, 2012.

[15] M. Rodrigues, K. Ezzedine, I. Hamzavi, A. G. Pandya, J. E. Harris, W. G. Vitiligo, "New discoveries in the pathogenesis and classification of vitiligo," Journal of the American Academy of Dermatology, vol. 77 no. 1,DOI: 10.1016/j.jaad.2016.10.048, 2017.

[16] I. C. Le Poole, A. Wankowicz-Kalinska, R. M. van den Wijngaard, B. J. Nickoloff, P. K. Das, "Autoimmune aspects of depigmentation in vitiligo," The Journal of Investigative Dermatology Symposium Proceedings, vol. 9 no. 1, pp. 68-72, DOI: 10.1111/j.1087-0024.2004.00825.x, 2004.

[17] Y. Jin, G. Andersen, D. Yorgov, T. M. Ferrara, S. Ben, K. M. Brownson, P. J. Holland, S. A. Birlea, J. Siebert, A. Hartmann, A. Lienert, N. van Geel, J. Lambert, R. M. Luiten, A. Wolkerstorfer, J. Wietze van der Veen, D. C. Bennett, A. Taïeb, K. Ezzedine, E. H. Kemp, D. J. Gawkrodger, A. P. Weetman, S. Kõks, E. Prans, K. Kingo, M. Karelson, M. R. Wallace, W. McCormack, A. Overbeck, S. Moretti, R. Colucci, M. Picardo, N. B. Silverberg, M. Olsson, Y. Valle, I. Korobko, M. Böhm, H. W. Lim, I. Hamzavi, L. Zhou, Q. S. Mi, P. R. Fain, S. A. Santorico, R. A. Spritz, "Genome-wide association studies of autoimmune vitiligo identify 23 new risk loci and highlight key pathways and regulatory variants," Nature Genetics, vol. 48 no. 11, pp. 1418-1424, DOI: 10.1038/ng.3680, 2016.

[18] A. Y. Lee, Y. H. Youm, N. H. Kim, H. Yang, W. I. Choi, "Keratinocytes in the depigmented epidermis of vitiligo are more vulnerable to trauma (suction) than keratinocytes in the normally pigmented epidermis, resulting in their apoptosis," The British Journal of Dermatology, vol. 151 no. 5, pp. 995-1003, DOI: 10.1111/j.1365-2133.2004.06136.x, 2004.

[19] T. W. Mak, W. C. Yeh, "Signaling for survival and apoptosis in the immune system," Arthritis Research, vol. 4, pp. S243-S252, DOI: 10.1186/ar569, 2002.

[20] A. Birol, U. Kisa, G. S. Kurtipek, F. Kara, M. Kocak, E. Erkek, O. Caglayan, "Increased tumor necrosis factor alpha (TNF- α ) and interleukin 1 alpha (IL1- α ) levels in the lesional skin of patients with nonsegmental vitiligo," International Journal of Dermatology, vol. 45 no. 8, pp. 992-993, DOI: 10.1111/j.1365-4632.2006.02744.x, 2006.

[21] Y. J. Suzuki, "Cell signaling pathways for the regulation of GATA4 transcription factor: implications for cell growth and apoptosis," Cellular Signalling, vol. 23 no. 7, pp. 1094-1099, DOI: 10.1016/j.cellsig.2011.02.007, 2011.

[22] S. Stromberg, M. G. Björklund, A. Asplund, R. Rimini, J. Lundeberg, P. Nilsson, F. Pontén, M. J. Olsson, "Transcriptional profiling of melanocytes from patients with vitiligo vulgaris," Pigment Cell & Melanoma Research, vol. 21 no. 2, pp. 162-171, DOI: 10.1111/j.1755-148X.2007.00429.x, 2008.

[23] N. Wang, H. Y. Tan, S. Li, Y. Feng, "Atg9b deficiency suppresses autophagy and potentiates endoplasmic reticulum stress-associated hepatocyte apoptosis in hepatocarcinogenesis," Theranostics, vol. 7 no. 8, pp. 2325-2338, DOI: 10.7150/thno.18225, 2017.

[24] C. Regazzetti, F. Joly, C. Marty, M. Rivier, B. Mehul, P. Reiniche, C. Mounier, Y. Rival, D. Piwnica, M. Cavalié, B. Chignon-Sicard, R. Ballotti, J. Voegel, T. Passeron, "Transcriptional analysis of vitiligo skin reveals the alteration of WNT pathway: a promising target for repigmenting vitiligo patients," The Journal of Investigative Dermatology, vol. 135 no. 12, pp. 3105-3114, DOI: 10.1038/jid.2015.335, 2015.

[25] R. Yu, R. Broady, Y. Huang, Y. Wang, J. Yu, M. Gao, M. Levings, S. Wei, S. Zhang, A. Xu, M. Su, J. Dutz, X. Zhang, Y. Zhou, "Transcriptome analysis reveals markers of aberrantly activated innate immunity in vitiligo lesional and non-lesional skin," PLoS One, vol. 7 no. 12, article e51040,DOI: 10.1371/journal.pone.0051040, 2012.

[26] Y. Yao, C. Dou, Z. Lu, X. Zheng, Q. Liu, "MACC1 suppresses cell apoptosis in hepatocellular carcinoma by targeting the HGF/c-MET/AKT pathway," Cellular Physiology and Biochemistry, vol. 35 no. 3, pp. 983-996, DOI: 10.1159/000369754, 2015.

[27] J. C. Samame Perez-Vargas, P. Biondani, C. Maggi, M. Gariboldi, A. Gloghini, A. Inno, C. Volpi, A. Gualeni, M. di Bartolomeo, F. de Braud, A. Castano, I. Bossi, F. Pietrantonio, "Role of cMET in the development and progression of colorectal cancer," International Journal of Molecular Sciences, vol. 14 no. 9, pp. 18056-18077, DOI: 10.3390/ijms140918056, 2013.

[28] D. A. Liebermann, B. Hoffman, "Gadd45 in the response of hematopoietic cells to genotoxic stress," Blood Cells, Molecules and Diseases, vol. 39 no. 3, pp. 329-335, DOI: 10.1016/j.bcmd.2007.06.006, 2007.

[29] L. Bouchier-Hayes, S. J. Martin, "CARD games in apoptosis and immunity," EMBO Reports, vol. 3 no. 7, pp. 616-621, DOI: 10.1093/embo-reports/kvf139, 2002.

[30] F. Shi, B. W. Kong, J. J. Song, J. Y. Lee, R. L. Dienglewicz, G. F. Erf, "Understanding mechanisms of vitiligo development in Smyth line of chickens by transcriptomic microarray analysis of evolving autoimmune lesions," BMC Immunology, vol. 13 no. 1, article 18,DOI: 10.1186/1471-2172-13-18, 2012.

[31] Y. Jin, S. A. Birlea, P. R. Fain, R. A. Spritz, "Genetic variations in NALP1 are associated with generalized vitiligo in a Romanian population," The Journal of Investigative Dermatology, vol. 127 no. 11, pp. 2558-2562, DOI: 10.1038/sj.jid.5700953, 2007.

[32] A. A. Beg, D. Baltimore, "An essential role for NF- κ B in preventing TNF- α -induced cell death," Science, vol. 274 no. 5288, pp. 782-784, DOI: 10.1126/science.274.5288.782, 1996.

[33] Y. L. Fan, L. Chen, J. Wang, Q. Yao, J. Q. Wan, "Over expression of PPP2R2C inhibits human glioma cells growth through the suppression of mTOR pathway," FEBS Letters, vol. 587 no. 24, pp. 3892-3897, DOI: 10.1016/j.febslet.2013.09.029, 2013.

[34] Y. Yamaguchi, A. Morita, A. Maeda, V. J. Hearing, "Regulation of skin pigmentation and thickness by Dickkopf 1 (DKK1)," The Journal of Investigative Dermatology Symposium Proceedings, vol. 14 no. 1, pp. 73-75, DOI: 10.1038/jidsymp.2009.4, 2009.

[35] J. Xie, D. B. Xiang, H. Wang, C. Zhao, J. Chen, F. Xiong, T. Y. Li, X. L. Wang, "Inhibition of Tcf-4 induces apoptosis and enhances chemosensitivity of colon cancer cells," PLoS One, vol. 7 no. 9, article e45617,DOI: 10.1371/journal.pone.0045617, 2012.

[36] R. K. Singh, K. M. Lee, I. Vujkovic-Cvijin, D. Ucmak, B. Farahnik, M. Abrouk, M. Nakamura, T. H. Zhu, T. Bhutani, M. Wei, W. Liao, "The role of IL-17 in vitiligo: a review," Autoimmunity Reviews, vol. 15 no. 4, pp. 397-404, DOI: 10.1016/j.autrev.2016.01.004, 2016.

[37] Y. Kotobuki, A. Tanemura, L. Yang, S. Itoi, M. Wataya-Kaneda, H. Murota, M. Fujimoto, S. Serada, T. Naka, I. Katayama, "Dysregulation of melanocyte function by Th17-related cytokines: significance of Th17 cell infiltration in autoimmune vitiligo vulgaris," Pigment Cell & Melanoma Research, vol. 25 no. 2, pp. 219-230, DOI: 10.1111/j.1755-148X.2011.00945.x, 2012.

[38] S. Singh, U. Singh, S. S. Pandey, "Serum concentration of IL-6, IL-2, TNF- α , and IFN γ in vitiligo patients," Indian Journal of Dermatology, vol. 57 no. 1, pp. 12-14, DOI: 10.4103/0019-5154.92668, 2012.

[39] H. S. Yu, K. L. Chang, C. L. Yu, H. F. Li, M. T. Wu, C. S. Wu, C. S. Wu, "Alterations in IL-6, IL-8, GM-CSF. TNF- α , and IFN- γ release by peripheral mononuclear cells in patients with active vitiligo," The Journal of Investigative Dermatology, vol. 108 no. 4, pp. 527-529, DOI: 10.1111/1523-1747.ep12289743, 1997.

[40] M. Singh, M. S. Mansuri, M. A. Parasrampuria, R. Begum, "Interleukin 1- α : a modulator of melanocyte homeostasis in vitiligo," Biochemistry & Analytical Biochemistry, vol. 5,DOI: 10.4172/2161-1009.1000273, 2016.

[41] K. R. Patel, V. Singam, S. Rastogi, H. H. Lee, N. B. Silverberg, J. I. Silverberg, "Association of vitiligo with hospitalization for mental health disorders in US adults," Journal of the European Academy of Dermatology and Venereology, vol. 33 no. 1, pp. 191-197, DOI: 10.1111/jdv.15255, 2019.

[42] P. M. Plonka, T. Passeron, M. Brenner, D. J. Tobin, S. Shibahara, A. Thomas, A. Slominski, A. L. Kadekaro, D. Hershkovitz, E. Peters, J. J. Nordlund, Z. Abdel-Malek, K. Takeda, R. Paus, J. P. Ortonne, V. J. Hearing, K. U. Schallreuter, "What are melanocytes really doing all day long…?," Experimental Dermatology, vol. 18 no. 9, pp. 799-819, DOI: 10.1111/j.1600-0625.2009.00912.x, 2009.

[43] S. Sarkar, T. Sarkar, A. Sarkar, S. Das, "Vitiligo and psychiatric morbidity: a profile from a vitiligo clinic of a rural-based tertiary care center of eastern India," Indian Journal of Dermatology, vol. 63 no. 4, pp. 281-284, DOI: 10.4103/ijd.IJD_142_18, 2018.

[44] P. Blanchet, M. Bebin, S. Bruet, G. M. Cooper, M. L. Thompson, B. Duban-Bedu, B. Gerard, A. Piton, S. Suckno, C. Deshpande, V. Clowes, J. Vogt, P. Turnpenny, M. P. Williamson, Y. Alembik, Clinical Sequencing Exploratory Research Study Consortium, Deciphering Developmental Disorders Consortium, E. Glasgow, A. McNeill, "MYT1L mutations cause intellectual disability and variable obesity by dysregulating gene expression and development of the neuroendocrine hypothalamus," PLoS Genetics, vol. 13 no. 8, article e1006957,DOI: 10.1371/journal.pgen.1006957, 2017.

[45] X. Liang, M. R. Song, Z. Xu, G. M. Lanuza, Y. Liu, T. Zhuang, Y. Chen, S. L. Pfaff, S. M. Evans, Y. Sun, "Isl1 is required for multiple aspects of motor neuron development," Molecular and Cellular Neurosciences, vol. 47 no. 3, pp. 215-222, DOI: 10.1016/j.mcn.2011.04.007, 2011.

[46] H. Zhang, W. P. Wang, T. Guo, J. C. Yang, P. Chen, K. T. Ma, Y. F. Guan, C. Y. Zhou, "The LIM-homeodomain protein ISL1 activates insulin gene promoter directly through synergy with BETA2," Journal of Molecular Biology, vol. 392 no. 3, pp. 566-577, DOI: 10.1016/j.jmb.2009.07.036, 2009.

[47] H. Shimomura, T. Sanke, T. Hanabusa, K. Tsunoda, H. Furuta, K. Nanjo, "Nonsense mutation of islet-1 gene (Q310X) found in a type 2 diabetic patient with a strong family history," Diabetes, vol. 49 no. 9, pp. 1597-1600, DOI: 10.2337/diabetes.49.9.1597, 2000.

[48] J. L. Munoz-Munoz, J. R. Acosta-Motos, F. Garcia-Molina, R. Varon, P. A. Garcia-Ruíz, J. Tudela, F. Garcia-Cánovas, J. N. Rodríguez-López, "Tyrosinase inactivation in its action on dopa," Biochimica et Biophysica Acta (BBA) - Proteins and Proteomics, vol. 1804 no. 7, pp. 1467-1475, DOI: 10.1016/j.bbapap.2010.02.015, 2010.

[49] A. M. Bolger, M. Lohse, B. Usadel, "Trimmomatic: a flexible trimmer for Illumina sequence data," Bioinformatics, vol. 30 no. 15, pp. 2114-2120, DOI: 10.1093/bioinformatics/btu170, 2014.

[50] E. Kopylova, L. Noe, H. Touzet, "SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data," Bioinformatics, vol. 28 no. 24, pp. 3211-3217, DOI: 10.1093/bioinformatics/bts611, 2012.

[51] I. Y. Zhbannikov, S. S. Hunter, J. A. Foster, M. L. Settles, "SeqyClean: a pipeline for high-throughput sequence data preprocessing," Proceedings of the 8th ACM International Conference on Bioinformatics, Computational Biology,and Health Informatics - ACM-BCB '17,DOI: 10.1145/3107411.3107446, .

[52] A. Dobin, C. A. Davis, F. Schlesinger, J. Drenkow, C. Zaleski, S. Jha, P. Batut, M. Chaisson, T. R. Gingeras, "STAR: ultrafast universal RNA-seq aligner," Bioinformatics, vol. 29 no. 1, pp. 15-21, DOI: 10.1093/bioinformatics/bts635, 2013.

[53] Y. Liao, G. K. Smyth, W. Shi, "featureCounts: an efficient general purpose program for assigning sequence reads to genomic features," Bioinformatics, vol. 30 no. 7, pp. 923-930, DOI: 10.1093/bioinformatics/btt656, 2014.

[54] M. E. Ritchie, B. Phipson, D. Wu, Y. Hu, C. W. Law, W. Shi, G. K. Smyth, "limma powers differential expression analyses for RNA-sequencing and microarray studies," Nucleic Acids Research, vol. 43 no. 7, article e47,DOI: 10.1093/nar/gkv007, 2015.

[55] D. W. Huang, B. T. Sherman, Q. Tan, J. R. Collins, W. G. Alvord, J. Roayaei, R. Stephens, M. W. Baseler, H. C. Lane, R. A. Lempicki, "The DAVID gene functional classification tool: a novel biological module-centric algorithm to functionally analyze large gene lists," Genome Biology, vol. 8 no. 9, article R183,DOI: 10.1186/gb-2007-8-9-r183, 2007.