Table of Contents Author Guidelines Submit a Manuscript
Oxidative Medicine and Cellular Longevity
Volume 2019, Article ID 2841814, 15 pages
Research Article

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

1Dermatology Research Unit, Division of Translational Medicine, Research Department, OPC 5th Floor, Sidra Medicine, Doha, PO Box 26999, Qatar
2Human Genetics, Translational Medicine, OPC 5th Floor, Sidra Medicine, Doha, PO Box 26999, Qatar

Correspondence should be addressed to Konduru Seetharama Sastry; gro.ardis@urudnoks and Aouatef Ismail Chouchane; gro.ardis@enahcuohca

Received 28 July 2019; Accepted 1 October 2019; Published 4 December 2019

Academic Editor: Demetrios Kouretas

Copyright © 2019 Konduru Seetharama Sastry et al. 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.


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.

1. Introduction

Reactive oxygen species (ROS) such as hydrogen peroxide (H2O2), 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 H2O2 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 H2O2-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 H2O2-induced oxidative stress in a time- and dose-dependent manner.

2. Results and Discussion

2.1. Generation of Oxidative Stress in Melanocytes Using H2O2 and Its Detection

Various batches of isolated primary human epidermal melanocytes (HEM) respond differently to same treatments with various drugs including H2O2. 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, H2O2 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 H2O2, cells start losing their dendriticity, a first hallmark of cell death process (Figure 1(d)).

Figure 1: Characterization of melanocytes and H2O2-induced oxidative death: (a) PIG1 and primary HEM were incubated with L-DOPA or PBS (control) for 5 h and images were obtained using phase contrast microscopy. Note that the melanocytes are branched and stained black by L-DOPA, confirming the presence of DOPA oxidase (tyrosinase) activity. (b) Indicated cells were cultured, and cell lysates were probed with MITF and TYRP-1 antibodies using western blotting. Equal loading was confirmed using β-actin antibodies. (c) PIG1 melanocytes were either left untreated or incubated with indicated concentration of H2O2 for 2 h at 37°C. The oxidative stress was detected by staining cells with CellROX Orange dye for 30 min, and images were obtained by the Evos fluorescent microscope. (d) To observe the morphological changes, phase contrast microscopic images were obtained 3 h after exposure to H2O2. (e) PIG1 melanocytes were cultured with indicated concentration of H2O2 for 6 h, 24 h, and 48 h. Cell viability was checked by trypan blue dye exclusion assay. Mean (three independent observations) values of viability percentages were plotted at different time intervals. Greater reduction in cell viability was found with increasing concentration of H2O2. Grey box includes the time points chosen for RNA-seq experiments. (f) Flow chart showing the different time points and concentrations of H2O2 used in RNA-seq experiments.
2.2. H2O2-Induced Cell Death in PIG1 Melanocytes

To choose the right concentration of H2O2 and dosing time points for RNA-seq experiments, we measured the cytotoxicity of H2O2 in melanocytes by exposing them to increasing concentrations of H2O2 for various exposure times. Cell viability was found to decrease in an inversely correlated manner to both dose and time of H2O2 exposure (Figure 1(e)). Across all concentrations of H2O2, 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 H2O2, it was reduced by 50% and 90% with 350 μM H2O2 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 H2O2 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 H2O2-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 H2O2 for 48 h. Similar pattern is observed with all other treatments (not shown).

Figure 2: (a) MDS plot: multidimensional scaling plot shows the clustering of controls and H2O2 treated samples. A large disparity in expression of genes between controls and H2O2 treated with various doses and time points was clearly evident. (b) Volcano plot: the transcriptome of cells treated with 250 μM H2O2 for 48 h was compared to control. The -axis represents the log2 fold changes (log2FC), and the -axis specifies the statistical significance (negative logarithm to the base 10 of FDR) of differential expression. Green vertical and horizontal lines reflect the filtering criteria ( and ). Red and blue dots represent up- and downregulated genes, respectively. (c) The Venn diagram shows the number of genes upregulated (Red) or downregulated (blue) by treatment with H2O2 for 6 h (a), 24 h (b), and 48 h (c) compared to controls. The numbers of differentially regulated genes in two or three groups are shown in intersections.
2.4. H2O2 Alters Gene Expression in a Time- and Dose-Dependent Manner in Melanocytes

H2O2 altered the expression changes in a time- and dose-dependent manner. We noticed a general trend that the number of genes induced by H2O2 exceeded the number of suppressed ones across time and concentrations (Figure 2(c)) (hypergeometric 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 H2O2 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.

Figure 3: (a) The hierarchical clustering of top 100 differentially expressed genes between 250 μM H2O2 at 48 h treatment and control. Red color is the upregulated and green color is the downregulated genes. (b) GO enrichment analysis of DEGs between 250 μM H2O2 samples and control samples at 48 h according to the DAVID bioinformatics tool. (c) Biological pathway enrichment analysis of DEGs between 250 μM H2O2 and control for 48 h using the DAVID bioinformatics tool.
Table 1: List of top 25 most upregulated genes by treatment with 250 μM H2O2 for 48 h. These genes passed the and abs.log2FC of 2.
Table 2: List of top 25 most downregulated genes by treatment with 250 μM H2O2 for 48 h. These genes passed the and abs.log2FC of 2.
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 H2O2-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 H2O2 for 48 h. These genes passed the and abs.log2FC of 2.
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 H2O2 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 H2O2-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.

Figure 4: (a) PIG1 cells and primary HEM were treated with 250 μM H2O2 for 48 h, and cell viability was measured by trypan blue dye exclusion assay. (b) PIG1 and primary HEM cells were treated as in (a), and mRNA extracted. Relative gene expression of indicated genes was measured using quantitative RT-PCR. GAPDH was used as internal control. (c) PIG1 cells and primary HEM were treated as in (a); cell lysates were collected and probed with indicated antibodies using western blotting. (d) The relative expression levels induced by treatment with 250 μM H2O2 for 48 h of BH family proapoptotic (blue) and antiapoptotic (green) are shown. Only those members whose value is <0.05 are shown.
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 [1216]. 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 H2O2-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 5: Flowchart depicts diverse signaling mechanisms altered by ROS in melanocytes. Red and green arrows indicate upregulation and downregulation of genes, respectively.

(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 H2O2 (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, H2O2 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 H2O2 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 H2O2. ATG9B, whose up expression was more prominent at 24 h after treatment with H2O2, is also well known to participate in autophagy [23]. Other autophagy genes either unchanged or downregulated may suggest that autophagy in H2O2-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 H2O2 (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 H2O2-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 H2O2 (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 H2O2 (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 H2O2 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 H2O2, 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 H2O2 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 H2O2 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 H2O2 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 H2O2. 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, H2O2 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 H2O2 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 H2O2 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 H2O2. 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 [4143], but no role of these TFs has been previously reported in melanocytes.

Table 4: Differentially regulated transcription factors between control and 250 μM H2O2 for 48 h and their putative functions. These genes passed the and abs.log2FC of 2.

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% CO2. 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 H2O2 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 log2 fold and . 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.


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.

Supplementary Materials

Supplementary Figure-1: comparison of the DEGs detected by three methods: Venn diagram of DEGs detected by the edgeR, DESeq, and limma. Supplementary Figure-2: the hierarchical clustering of top 100 differentially expressed genes (A): control vs. 100 μM H2O2 (B): control vs. 350 μM H2O2 at 48 h treatment and control. Red color is the upregulated and green color is the downregulated genes. Supplementary Figure-3: GO enrichment analysis of DEGs (A): 100 μM H2O2 and control samples at 48 h; (B): 350 μM H2O2 and control samples at 48 h, according to the DAVID bioinformatics tool. Supplementary Figure-4: biological pathway enrichment analysis of DEGs (A): 100 μM H2O2 and control for 48 h; (B): 350 μM H2O2 and control for 48 h using the DAVID bioinformatics tool. Supplementary Figure-5: PIG1 cells were treated with 250 μM H2O2 for 48 h, and cell viability was measured by trypan blue dye exclusion assay (A). Relative gene expression of indicated genes was measured using quantitative RT-PCR (B). GAPDH and beta-actin were used as internal controls separately. Supplementary Table-1: list of differentially regulated selected apoptotic proteins after 24 h of treatment with H2O2. (Supplementary Materials)


  1. M. Schieber and N. S. Chandel, “ROS function in redox signaling and oxidative stress,” Current Biology, vol. 24, no. 10, pp. R453–R462, 2014. View at Publisher · View at Google Scholar · View at Scopus
  2. L. Denat, A. L. Kadekaro, L. Marrot, S. A. Leachman, and Z. A. Abdel-Malek, “Melanocytes as instigators and victims of oxidative stress,” The Journal of Investigative Dermatology, vol. 134, no. 6, pp. 1512–1518, 2014. View at Publisher · View at Google Scholar · View at Scopus
  3. A. Jain, J. Mal, V. Mehndiratta, R. Chander, and S. K. Patra, “Study of oxidative stress in vitiligo,” Indian Journal of Clinical Biochemistry, vol. 26, no. 1, pp. 78–81, 2011. View at Publisher · View at Google Scholar · View at Scopus
  4. K. U. Schallreuter, J. Moore, J. M. Wood et al., “In vivo and in vitro evidence for hydrogen peroxide (H2O2) 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. View at Google Scholar
  5. K. Jimbow, H. Chen, J. S. Park, and 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, 2001. View at Publisher · View at Google Scholar · View at Scopus
  6. O. A. Arowojolu, S. J. Orlow, N. Elbuluk, and 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, 2017. View at Publisher · View at Google Scholar · View at Scopus
  7. J. A. Mosenson, K. Flood, J. Klarquist et al., “Preferential secretion of inducible HSP70 by vitiligo melanocytes under stress,” Pigment Cell & Melanoma Research, vol. 27, no. 2, pp. 209–220, 2014. View at Publisher · View at Google Scholar · View at Scopus
  8. G. Yang, G. Zhang, M. R. Pittelkow, M. Ramoni, and 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, 2006. View at Publisher · View at Google Scholar · View at Scopus
  9. I. C. Le Poole, F. M. van den Berg, R. M. van den Wijngaard et al., “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, 1997. View at Publisher · View at Google Scholar · View at Scopus
  10. Z. Zhou, C. Y. Li, K. Li, T. Wang, B. Zhang, and 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, 2009. View at Publisher · View at Google Scholar · View at Scopus
  11. P. Wang, Y. Li, H. Nie et al., “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, 2016. View at Publisher · View at Google Scholar · View at Scopus
  12. A. M. Abdel-Aal, M. A. Kasem, and 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, pp. 1–11, 2002. View at Google Scholar
  13. C. L. Huang, J. J. Nordlund, and R. Boissy, “Vitiligo: a manifestation of apoptosis?” American Journal of Clinical Dermatology, vol. 3, no. 5, pp. 301–308, 2002. View at Publisher · View at Google Scholar · View at Scopus
  14. R. Kumar and D. Parsad, “Melanocytorrhagy and apoptosis in vitiligo: connecting jigsaw pieces,” Indian Journal of Dermatology, Venereology and Leprology, vol. 78, no. 1, pp. 19–23, 2012. View at Publisher · View at Google Scholar · View at Scopus
  15. M. Rodrigues, K. Ezzedine, I. Hamzavi, A. G. Pandya, J. E. Harris, and W. G. Vitiligo, “New discoveries in the pathogenesis and classification of vitiligo,” Journal of the American Academy of Dermatology, vol. 77, no. 1, pp. 1–13, 2017. View at Publisher · View at Google Scholar · View at Scopus
  16. I. C. Le Poole, A. Wankowicz-Kalinska, R. M. van den Wijngaard, B. J. Nickoloff, and P. K. Das, “Autoimmune aspects of depigmentation in vitiligo,” The Journal of Investigative Dermatology Symposium Proceedings, vol. 9, no. 1, pp. 68–72, 2004. View at Publisher · View at Google Scholar · View at Scopus
  17. Y. Jin, G. Andersen, D. Yorgov et al., “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, 2016. View at Publisher · View at Google Scholar · View at Scopus
  18. A. Y. Lee, Y. H. Youm, N. H. Kim, H. Yang, and 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, 2004. View at Publisher · View at Google Scholar · View at Scopus
  19. T. W. Mak and W. C. Yeh, “Signaling for survival and apoptosis in the immune system,” Arthritis Research, vol. 4, Supplement 3, pp. S243–S252, 2002. View at Publisher · View at Google Scholar
  20. A. Birol, U. Kisa, G. S. Kurtipek et al., “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, 2006. View at Publisher · View at Google Scholar · View at Scopus
  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, 2011. View at Publisher · View at Google Scholar · View at Scopus
  22. S. Stromberg, M. G. Björklund, A. Asplund et al., “Transcriptional profiling of melanocytes from patients with vitiligo vulgaris,” Pigment Cell & Melanoma Research, vol. 21, no. 2, pp. 162–171, 2008. View at Publisher · View at Google Scholar · View at Scopus
  23. N. Wang, H. Y. Tan, S. Li, and Y. Feng, “Atg9b deficiency suppresses autophagy and potentiates endoplasmic reticulum stress-associated hepatocyte apoptosis in hepatocarcinogenesis,” Theranostics, vol. 7, no. 8, pp. 2325–2338, 2017. View at Publisher · View at Google Scholar · View at Scopus
  24. C. Regazzetti, F. Joly, C. Marty et al., “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, 2015. View at Publisher · View at Google Scholar · View at Scopus
  25. R. Yu, R. Broady, Y. Huang et al., “Transcriptome analysis reveals markers of aberrantly activated innate immunity in vitiligo lesional and non-lesional skin,” PLoS One, vol. 7, no. 12, article e51040, 2012. View at Publisher · View at Google Scholar · View at Scopus
  26. Y. Yao, C. Dou, Z. Lu, X. Zheng, and 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, 2015. View at Publisher · View at Google Scholar · View at Scopus
  27. J. C. Samame Perez-Vargas, P. Biondani, C. Maggi et al., “Role of cMET in the development and progression of colorectal cancer,” International Journal of Molecular Sciences, vol. 14, no. 9, pp. 18056–18077, 2013. View at Publisher · View at Google Scholar · View at Scopus
  28. D. A. Liebermann and B. Hoffman, “Gadd45 in the response of hematopoietic cells to genotoxic stress,” Blood Cells, Molecules and Diseases, vol. 39, no. 3, pp. 329–335, 2007. View at Publisher · View at Google Scholar · View at Scopus
  29. L. Bouchier-Hayes and S. J. Martin, “CARD games in apoptosis and immunity,” EMBO Reports, vol. 3, no. 7, pp. 616–621, 2002. View at Publisher · View at Google Scholar · View at Scopus
  30. F. Shi, B. W. Kong, J. J. Song, J. Y. Lee, R. L. Dienglewicz, and 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, 2012. View at Publisher · View at Google Scholar · View at Scopus
  31. Y. Jin, S. A. Birlea, P. R. Fain, and 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, 2007. View at Publisher · View at Google Scholar · View at Scopus
  32. A. A. Beg and D. Baltimore, “An essential role for NF-κB in preventing TNF-α-induced cell death,” Science, vol. 274, no. 5288, pp. 782–784, 1996. View at Publisher · View at Google Scholar · View at Scopus
  33. Y. L. Fan, L. Chen, J. Wang, Q. Yao, and 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, 2013. View at Publisher · View at Google Scholar · View at Scopus
  34. Y. Yamaguchi, A. Morita, A. Maeda, and 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, 2009. View at Publisher · View at Google Scholar · View at Scopus
  35. J. Xie, D. B. Xiang, H. Wang et al., “Inhibition of Tcf-4 induces apoptosis and enhances chemosensitivity of colon cancer cells,” PLoS One, vol. 7, no. 9, article e45617, 2012. View at Publisher · View at Google Scholar · View at Scopus
  36. R. K. Singh, K. M. Lee, I. Vujkovic-Cvijin et al., “The role of IL-17 in vitiligo: a review,” Autoimmunity Reviews, vol. 15, no. 4, pp. 397–404, 2016. View at Publisher · View at Google Scholar · View at Scopus
  37. Y. Kotobuki, A. Tanemura, L. Yang et al., “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, 2012. View at Publisher · View at Google Scholar · View at Scopus
  38. S. Singh, U. Singh, and 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, 2012. View at Publisher · View at Google Scholar · View at Scopus
  39. H. S. Yu, K. L. Chang, C. L. Yu et al., “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, 1997. View at Publisher · View at Google Scholar
  40. M. Singh, M. S. Mansuri, M. A. Parasrampuria, and R. Begum, “Interleukin 1-α: a modulator of melanocyte homeostasis in vitiligo,” Biochemistry & Analytical Biochemistry, vol. 5, 2016. View at Publisher · View at Google Scholar
  41. K. R. Patel, V. Singam, S. Rastogi, H. H. Lee, N. B. Silverberg, and 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, 2019. View at Publisher · View at Google Scholar · View at Scopus
  42. P. M. Plonka, T. Passeron, M. Brenner et al., “What are melanocytes really doing all day long…?” Experimental Dermatology, vol. 18, no. 9, pp. 799–819, 2009. View at Publisher · View at Google Scholar · View at Scopus
  43. S. Sarkar, T. Sarkar, A. Sarkar, and 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, 2018. View at Publisher · View at Google Scholar · View at Scopus
  44. P. Blanchet, M. Bebin, S. Bruet et al., “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, 2017. View at Publisher · View at Google Scholar · View at Scopus
  45. X. Liang, M. R. Song, Z. Xu et al., “Isl1 is required for multiple aspects of motor neuron development,” Molecular and Cellular Neurosciences, vol. 47, no. 3, pp. 215–222, 2011. View at Publisher · View at Google Scholar · View at Scopus
  46. H. Zhang, W. P. Wang, T. Guo et al., “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, 2009. View at Publisher · View at Google Scholar · View at Scopus
  47. H. Shimomura, T. Sanke, T. Hanabusa, K. Tsunoda, H. Furuta, and 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, 2000. View at Publisher · View at Google Scholar · View at Scopus
  48. J. L. Munoz-Munoz, J. R. Acosta-Motos, F. Garcia-Molina et al., “Tyrosinase inactivation in its action on dopa,” Biochimica et Biophysica Acta (BBA) - Proteins and Proteomics, vol. 1804, no. 7, pp. 1467–1475, 2010. View at Publisher · View at Google Scholar · View at Scopus
  49. A. M. Bolger, M. Lohse, and B. Usadel, “Trimmomatic: a flexible trimmer for Illumina sequence data,” Bioinformatics, vol. 30, no. 15, pp. 2114–2120, 2014. View at Publisher · View at Google Scholar · View at Scopus
  50. E. Kopylova, L. Noe, and H. Touzet, “SortMeRNA: fast and accurate filtering of ribosomal RNAs in metatranscriptomic data,” Bioinformatics, vol. 28, no. 24, pp. 3211–3217, 2012. View at Publisher · View at Google Scholar · View at Scopus
  51. I. Y. Zhbannikov, S. S. Hunter, J. A. Foster, and M. L. Settles, “SeqyClean: a pipeline for high-throughput sequence data preprocessing,” in Proceedings of the 8th ACM International Conference on Bioinformatics, Computational Biology,and Health Informatics - ACM-BCB '17, Boston, Massachusetts, USA, 2017. View at Publisher · View at Google Scholar · View at Scopus
  52. A. Dobin, C. A. Davis, F. Schlesinger et al., “STAR: ultrafast universal RNA-seq aligner,” Bioinformatics, vol. 29, no. 1, pp. 15–21, 2013. View at Publisher · View at Google Scholar · View at Scopus
  53. Y. Liao, G. K. Smyth, and W. Shi, “featureCounts: an efficient general purpose program for assigning sequence reads to genomic features,” Bioinformatics, vol. 30, no. 7, pp. 923–930, 2014. View at Publisher · View at Google Scholar · View at Scopus
  54. M. E. Ritchie, B. Phipson, D. Wu et al., “limma powers differential expression analyses for RNA-sequencing and microarray studies,” Nucleic Acids Research, vol. 43, no. 7, article e47, 2015. View at Publisher · View at Google Scholar · View at Scopus
  55. D. W. Huang, B. T. Sherman, Q. Tan et al., “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, 2007. View at Publisher · View at Google Scholar · View at Scopus