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

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


Introduction
Reactive oxygen species (ROS) such as hydrogen peroxide (H 2 O 2 ), 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 2 O 2 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 2 O 2induced transcriptional changes in melanocytes are known, a comprehensive view of stress-induced gene regulation is still elusive. While much of the understanding of ROSinduced 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 2 O 2 -induced oxidative stress in a time-and dose-dependent manner.

Generation of Oxidative Stress in Melanocytes
Using H 2 O 2 and Its Detection. Various batches of isolated primary human epidermal melanocytes (HEM) respond differently to same treatments with various drugs including H 2 O 2 . 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 2 O 2 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 2 O 2 , cells start losing their dendriticity, a first hallmark of cell death process (Figure 1(d)).  (Figure 1(f)).

Gene Expression
Pattern. MDS plot shows a clear clustering of controls and H 2 O 2 -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 2 O 2 for 48 h. Similar pattern is observed with all other treatments (not shown).  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.

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 2 O 2 -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.
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 2 O 2 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 2 O 2 -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 medi-ated by the differential molecular abnormalities of these proteins.  vitiligo [12][13][14][15][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.

Cell Death
(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 RNAseq results suggest that the regulation of apoptosis in stressinduced 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 2 O 2 -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.
(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 2 O 2 ( 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 2 O 2 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 2 O 2 as seen in our study, an observation consistent with a previous report showing the downregulation of GATA3 in vitiligo melanocytes [22], may likely be 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 2 O 2 -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  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 2 O 2 (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 2 O 2 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 2 O 2 , 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.    receptors in melanocytes treated with H 2 O 2 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 2 O 2 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 2 O 2 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 2 O 2 . 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 2 O 2 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 2 O 2 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 2 O 2 include ZEB1, SNAI1, MYT1L, SOX21, KLF4, and GRHL2, all of which are well known to control epithelial to mesenchymal transition, development, or pluripotency.

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 2 O 2 . 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][42][43], but no role of these TFs has been previously reported in melanocytes.
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 poten-tially 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 doseand 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, failproof, 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.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 2 . HaCaT cells were cultured in M154 medium supplemented with HKGS. L-DOPA staining of melanocytes was performed as described previously [48].

Materials and Methods
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 2 O 2 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 2 fold change ðlog 2 FCÞ > 2:0 and 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.

Data Availability
The RNA-seq data used to support the findings of this study are available from the corresponding author upon request.

Conflicts of Interest
All authors declare no conflicts of interests.