Differential and Overlapping Effects of Melatonin and Its Metabolites on Keratinocyte Function: Bioinformatics and Metabolic Analyses

We investigated the effects of melatonin and its selected metabolites, i.e., N1-Acetyl-N2-formyl-5-methoxykynurenamine (AFMK) and 6-hydroxymelatonin (6(OH)Mel), on cultured human epidermal keratinocytes (HEKs) to assess their homeostatic activities with potential therapeutic implications. RNAseq analysis revealed a significant number of genes with distinct and overlapping patterns, resulting in common regulation of top diseases and disorders. Gene Set Enrichment Analysis (GSEA), Reactome FIViZ, and Ingenuity Pathway Analysis (IPA) showed overrepresentation of the p53-dependent G1 DNA damage response gene set, activation of p53 signaling, and NRF2-mediated antioxidative pathways. Additionally, GSEA exhibited an overrepresentation of circadian clock and antiaging signaling gene sets by melatonin derivatives and upregulation of extension of telomere signaling in HEKs, which was subsequently confirmed by increased telomerase activity in keratinocytes, indicating possible antiaging properties of metabolites of melatonin. Furthermore, Gene Ontology (GO) showed the activation of a keratinocyte differentiation program by melatonin, and GSEA indicated antitumor and antilipidemic potential of melatonin and its metabolites. IPA also indicated the role of Protein Kinase R (PKR) in interferon induction and antiviral response. In addition, the test compounds decreased lactate dehydrogenase A (LDHA) and lactate dehydrogenase C (LDHC) gene expression. These results were validated by qPCR and by Seahorse metabolic assay with significantly decreased glycolysis and lactate production under influence of AFMK or 6(OH)Mel in cells with a low oxygen consumption rate. In summary, melatonin and its metabolites affect keratinocytes’ functions via signaling pathways that overlap for each tested molecule with some distinctions.

To provide insight into the mechanisms of the above antioxidative, protective, metabolic, prodifferentiation, antiaging, anti-inflammatory, and antitumor mechanisms in the epidermis, we conducted RNAseq analyses using human epidermal keratinocytes treated with melatonin and its indolic (6-hydroxymelatonin) and kynuric N 1 -Acetyl-N 2 -formyl-5-methoxykynurenamine (AFMK) metabolites. We found overlapping and differential effects that were further investigated using biochemical and molecular assays. By using the comprehensive Omics method of estimation of gene expression, we obtained an insight not only into the expression of particular genes involved in the above processes, but also the correlation between these genes and their functional clustering. Based on these analyses, the ability of melatonin's metabolites to regulate different pathways affecting the physiological and pathological states of epidermis remains under discussion.  Test media was supplemented with 2 mM Glu-taMAX™ Supplement (L-alanyl-L-glutamine dipeptide in NaCl), Proteome Profiler Array kit (R&D Systems, Minneapolis, MN, USA), lactate colorimetric assay kit (Cell Biolabs, Inc. San Diego, CA, USA), TeloTAGGG Telomerase Elisa assay (Roche, Basel, Switzerland), and a set of human primers for real-time PCR (Eurofins Genomics, Ebersberg, Germany).

Cell Culture and Treatment
Neonatal Human Epidermal Keratinocytes (HEKn, primary cells) were isolated from foreskins collected at the Woman and Infant Hospital (UAB) and cultured on Petri dishes (TPP) in EpiGRO™ Human Epidermal Keratinocyte Complete Culture Media Kit (Millipore Merck KGaA, Darmstadt, Germany) as described previously [12]. Cells were seeded into 24-well plates, maintained until they reached 70% confluency, and stimulated for 24 h with melatonin, 6(OH)Mel or AFMK at the final concentration of 10 −5 M or with solvent only (control, 0.1% EtOH). HaCaT cells were grown in Dulbecco's Modified Eagle Medium (DMEM) media containing 5% charcoal stripped FBS followed by treatment with 10 −4 or 10 −5 M melatonin and 10 −5 M 6(OH)Mel and 2-hydroxymelatonin (2(OH)Mel) at 37 • C under 5% CO 2 for 6 and 24 h. We have used similar doses of ligands in our previous studies, which demonstrated the paracrine mechanism of action secondary to production and metabolism of melatonin in the epidermis [12,13,29,30]. After 6 and 24 h, cells were collected in lysis reagent and RNA isolated using Absolutely RNA Miniprep Kit (Agilent Technologies, Santa Clara, CA, USA).

Reverse Transcription Reaction and Real-Time PCR
To study gene expression in material acquired from HaCaT cells, cDNA was synthesized using 0.

RNA Sequencing and Bioinformatics Analysis
At least 200 ng of RNA from each HEKn sample with purity of 2.0 (A 260/280 ) was sent for RNA sequencing by Novogene (Sacramento, CA, USA). Prior to sending, the sample quality and concentration were determined using Cytation 5 Imaging Reader. After sequencing, the raw sequence FASTQ files were trimmed to remove primer adapters using Trim Galore version 0.6.2 (parameters used: -paired; -trim1). The trimmed FASTQ files were then aligned to the Gencode GRCh38 p7 Release 25 genome using STAR version 2.7.1a (parameters used: -outReadsUnmappedFastx; -outSAMtype BAM SortedByCoordinate; -outSAMattributes All). Transcript abundances were then calculated from the alignments using HTSeq count version 0.11.2 (parameters used: -m union; -r pos; -t exon; -igene_id; -a 10; -s no; -f bam). The raw counts were then loaded into DESeq2 using their default settings to normalize and perform pairwise differential expression. All raw data were deposited in NCBI GEO (GSE147588).

Bioinformatics
Pathway analysis was performed using Reactome and Ingenuity Pathway Analysis software (Ingenuity ® Systems). Gene set enrichment analysis (GSEA) and diagram analysis with overlapping gene scores into pathways were made using Reactome FIViZ software. Data analyzed using IPA software were processed with the following steps: for generating networks, a data set containing gene identifiers and corresponding expression values were uploaded into the application. Each identifier was mapped to its corresponding object in Ingenuity's Knowledge Base. A fold change cutoff of +/−2 was set to identify molecules whose expressions were significantly differentially regulated. These molecules, called Network Eligible molecules, were overlaid onto a global molecular network developed from information contained in Ingenuity's Knowledge Base. Networks of Network Eligible Molecules were then algorithmically generated based on their connectivity. The functional analysis identified the biological functions and/or diseases that were most significant to the entire data set. Molecules from the data set that met the fold change cutoff of +/−2 and were associated with biological functions and/or diseases in Ingenuity's Knowledge Base were considered for the analysis. Right-tailed Fisher's exact test was used to calculate the p value determining the probability that each biological function and/or disease assigned to that data set is due to chance alone.

Analysis of Metabolic Function
Primary human keratinocytes were isolated from skin biopsies of healthy controls as previously described [44]. Keratinocytes were seeded at 0.25 × 10 5 cells/well overnight on Seahorse 96-well XFe96 microplates in culture medium. The next day, cells were treated with melatonin, AFMK and 6(OH)Mel for 24 h, all at a final concentration of 10 −5 M. Just before the start of the experiment, cells were depleted of glucose for 1 h in a 37 • C non-CO 2 incubator. Oxygen consumption rate (OCR or mitochondrial respiration) and extracellular acidification rate (ECAR or glycolytic function) in live cells in real-time were measured using the XF Cell Mito Stress Test and Agilent Seahorse XF96 Extracellular Flux Analyzer (Seahorse Bioscience, North Billerica, MA, USA). The Cell Mito Stress Test media was supplemented with 2 mM GlutaMAX™ Supplement, 1 mM sodium pyruvate, and 25 mM glucose and injected according to the Mito Stress Test. Results are automatically generated from Wave data that was exported to Excel.

Proteome Profiler and L-Lactate Functional Assays
To measure cytokine release, conditioned media from HEKn cells collected after treatment with melatonin, AFMK or 6(OH)Mel for 24 h were submitted for human cytokine array assay using Proteome Profiler Array kit (R&D Systems). To evaluate L-lactate production in HEKn cell media, the Lactate Assay kit (Colorimetric, Cell Biolabs) was used and the values were recorded using Cytation 5 Imaging Reader. Telomerase activity was measured using the TeloTAGGG Telomerase Elisa assay (Roche) in the cell pellets obtained after 24 h treatment of HaCaT keratinocytes with AFMK, melatonin or 6(OH)Mel at a concentration of 10 −4 or 10 −5 M. The data were collected with Cytation 5 Imaging Reader. All assays followed procedures provided by manufacturers.

Statistical Analysis
Experiments were performed at least three times, with results expressed in each case as the mean + standard deviation (SD). Significant differences between results were determined by both by Mann-Whitney U test and Student t-test, and an appropriate post hoc analysis using GraphPad Prism 7.05 software (La Jolla, CA, USA). p value of less than 0.05 was considered statistically significant.

Bioinformatics Analysis
IPA demonstrated common regulation of top three diseases and disorders including cancer, dermatological diseases and conditions and organismal injuries and abnormalities by all compounds, with two additional ones having different sequence orders: reproductive system disease and endocrine system disorders for melatonin, reproductive system disease and psychological disorders for 6(OH)Mel, and endocrine system disorders and gastrointestinal diseases for AFMK (Table 1, Supplementary Materials Table S1). There were also similarities with some differences in regulation of molecular and cellular functions that included cell-to-cell signaling and interactions, gene expression, cell signaling, molecular transport and vitamin and mineral metabolisms for melatonin; molecular transport, cellto-cell signaling and interaction, cell signaling, vitamin and mineral metabolisms and cell morphology for 6(OH)Mel; molecular transport, cellular function and maintenance, cell signaling, vitamin and mineral metabolisms and cellular movement for AFMK (Table 1).
GSEA and Reactome FIViZ revealed upregulation of gene sets connected with responses to oxidative stress and DNA repair ( Table 2). IPA was consistent with this analysis by showing stimulation of the p53 pathway and p38 MAPK signaling, NER pathway, NRF2-mediated oxidative stress response by melatonin, AFMK and 6(OH)Mel and the role of BRCA1 in DNA damage response only by melatonin and AFMK ( Figure 2; Figure 3).
The profile and Z-scores of expressions of the genes involved in NRF2 downstream signaling affected by tested compounds are shown in Figure 3A. Expression of solute carrier family 6 member 9 (SLC6A9) was stimulated by melatonin treatment, inhibited by 6(OH)Mel and did not change after AFMK treatment. Expression of early growth response protein 1 (EGR1) was upregulated by melatonin, downregulated by 6(OH)Mel and showed no change for AFMK, carboxylesterase 2-phase 1 protein (CES2), EPH receptor A2 (EPHA2, which belongs to ephrin receptor subfamily), biliverdin reductase B (BLVRB), and transforming growth factor α (TGFA) and had similar expression levels to the control after treatment with AFMK or melatonin and was downregulated by 6(OH)Mel. SQSTM1 (sequestosome 1) had a similar Z-score for AFMK and 6(OH)Mel and the control and a negative Z-score value (-1.5) for melatonin. Both heparin binding EGF-like growth factor (HBEGF) and heme oxygenase 1 (HMOX1) genes were inhibited by the tested compounds. Expression of solute carrier family 7 member 11 (SLC7A11) and solute carrier family 2 member 3 (SLC2A3) was inhibited by melatonin and 6(OH)Mel, but not by AFMK. Thioredoxin reductase 1 (TXNRD1) was strongly upregulated by AFMK and to a lower degree by 6(OH)Mel and slightly by melatonin. This observation is consistent with reported role of AFMK in protection against oxidative stress [10,26,44]. Expression of neuregulin 1 (NRG1) was only stimulated by AFMK. Transporter solute carrier family 39 member 7 (SLC39A7) only had a positive Z-score for AFMK, and solute carrier family 39 member 4 (SLC39A4) only had a positive Z-score for 6(OH)Mel and melatonin. MAF bZIP transcription factor F (MAFF) and v-maf avian musculoaponeurotic fibrosarcoma oncogene homolog G (MAFG), crucial regulators of mammalian gene expression, were strongly upregulated only by melatonin with weak stimulation by AFMK and 6(OH)Mel for MAFG only. In summary, these analyses are consistent with previous experimental data showing that melatonin and its metabolites differentially regulate the Nrf2 signaling pathway in keratinocytes [12,23,27,36] and such regulation would depend on melatonin metabolism [22].   The profile and Z-scores of expressions of the genes involved in TP53 signaling are presented in Figure 3B. TP53BP1 and TP53TG1 were upregulated by AFMK. The expression levels of growth arrest and DNA damage-inducible protein gamma (GADD45G) and GADD45A increased after melatonin treatment, and TP53I3 and TP53INP2 were upregulated by AFMK and 6(OH)Mel. The increased GADD45G level is consistent with melatonin involvement in cellular stress responses and DNA repair processes, as well as tumor suppression [1,10,12,21,26,27].
Increased GADD45A expression by melatonin was associated with DNA damageinduced transcription of this gene mediated by both p53-dependent and p53-independent mechanisms and responses to environmental stress by mediating activation of the p38/JNK pathway via MTK1/MEKK4 kinase. P53I11 and TP53INP1 were also upregulated by AFMK and 6(OH)Mel. All compounds increased the expression of GADD45B, indicating their involvement in the response to environmental stress by activation of p38/JNK signaling [45].     The expression profile of genes involved in a mitochondrial metabolism was similar for all tested compounds (Table 3). However, minor differences were seen in the regulation of pathway of glycogen synthesis, which was upregulated by melatonin and AFMK, but not by 6(OH)Mel. The mitochondrial fatty acid beta-oxidation signaling was overrepresented only after melatonin treatment (Table 3) Table 4. GSEA shows an increase in extension of telomeres, telomere maintenance, telomere C-strand (lagging strand) synthesis and packaging of telomere end signaling, which indicates antiaging properties of melatonin, AFMK and 6(OH)Mel in the skin. These analyses are consistent with the recognized antiaging properties of melatonin [46], supported further by the functional data presented in Section 3.4. In addition, the GSEA also suggests an antiviral activity of melatonin and its metabolites in keratinocytes with upregulation of innate immunity. This is consistent with IPA, indicating a role of RIG1-like receptors in the antiviral innate immunity pathway for 6(OH)Mel as well as the role of PKR in the interferon induction and antiviral responses. The role of RIG1-like receptors in antiviral innate immunity and the role of PKR in interferon induction and antiviral response signaling were also seen after AFMK treatment. The same set of pathways was found in IPA for melatonin. These data are consistent with previous reports on antiviral properties of melatonin [38]. Table 4. Therapeutic implications for melatonin, AFMK and 6(OH)Mel based on GSEA. NES-normalized enriched score; FDR-false discovery rate; (×)-the effect is absent.
We also found an influence of melatonin on an expression of genes involved in the regulation of molecular clock in the skin through GSEA. We noticed an overrepresentation of circadian clock genes set. Differential expression analysis of RNAseq data indicated changes in expression of particular genes associated with circadian clock under influence of melatonin. There were upregulations of the NPAS2 (fold change = 2.3) and PPARγ (fold change = 2.7) genes as well as of the HDAC4 gene (fold change = 2.1) ( Figure 1A). The activation of PPARγ leads to strong anti-inflammatory effects [48] −1.1). Therefore, we suggest that melatonin can also affect the circadian molecular clock in keratinocytes in a similar manner as in other organs [49].
The epidermal lipids produced in keratinocytes play an essential role in the skin's barrier formation and function [50]. These lipids constitute a barrier against the loss of water and electrolytes and a barrier against microbial and viral invasion. However, elevated sebum excretion with lipid content is a key factor involved in the acne pathology [50,51]. The GSEA also suggested hypolipidemic and antiatherogenic potential for melatonin and its metabolites (Table 5). Table 5. Implications for cardiovascular system based on the GSEA of gene expression pattern affected by melatonin, AFMK and 6(OH)Mel in primary human epidermal keratinocytes. NES-normalized enriched score; (×)-the effect is absent. Interestingly, melatonin and 6(OH)Mel increased gamma-aminobutyric acid (GABA) receptor signaling in canonical IPA (Figure 2). The GABA receptors respond to GABA, the main inhibitory neurotransmitter in the mature vertebrate central nervous system (CNS). Therefore, we speculate that melatonin, AFMK and 6(OH)Mel might act similarly to GABA receptor agonists, which are the class of drugs that typically possesses sedative action with anticonvulsant, anxiolytic and muscle relaxant activities [52].

Reactome
Finally, the analysis indicates the anticancer potential of melatonin, AFMK and 6(OH)Mel by inhibition of pathways involved in process of carcinogenesis or by activation of anticancer immune signaling (Table 6). Antitumor properties of the tested compound were confirmed by an additional pathway analysis of differential expressed genes using the reactome platform. This showed negative regulation of Notch 4 signaling (with entities ratio = 0.008, p value = 0.04, FDR = 0.451; reactions ratio = 4.92). These analyses indicate complex action of melatonin and its metabolites against multiple types of tumors and are consistent with well-established anticancer activity of melatonin [1,8]. Interest-ingly, 6(OH)Mel and AFMK show inhibition of the osteoarthritis pathway in IPA with a Z-score = −0.535 and ratio = 0.102 and Z-score = −0.688 and ratio = 0.124, respectively. This analysis is consistent with the reported protective role caused by application of melatonin in experimental models of osteoarthritis [53][54][55][56]. Table 6. Indications for protective effects of melatonin, AFMK and 6(OH)Mel against common human tumors based on Canonical Ingenuity Pathway Analysis (IPA) of data obtained from primary human epidermal keratinocytes. (×)-the effect is absent.

Cytokine Array Analysis
Keratinocytes are the main source of cytokines in the epidermis and secrete CSF, TNFα, IL-1, IL-6, IL-3, IL-8, TGFα, TGFβ, and PDGF [57]. They can act as a complex regulatory network in the epidermis influenced by physiological or pathological signals [58,59]. The cytokines' activity is mediated by interaction with corresponding high affinity receptors [57]. The lower level of the expression of macrophage migration inhibitory factor (MIF) was noticed in DESeq2 analysis of RNAseq data for all tested compounds (fold change: 6(OH)Mel = −1.08, melatonin = −1.06, AFMK = −1.14). Tests using the cytokine profiler array have confirmed these results ( Figure 4A). We also noticed elevated interleukin 1 receptor antagonist (IL1ra) protein expression after treatment with melatonin and with 6(OH)Mel ( Figure 4B).
IL1ra is a competitive inhibitor of IL1, which participates in the process of silencing of inflammatory responses [60]. This suggests that melatonin and its derivatives can attenuate the proinflammatory conditions by increasing expression of ILra protein. Serpine family E member 1 (SERPINE1)'s protein expression was significantly elevated by 6(OH)Mel only ( Figure 4C). SERPINE1 is a protein required for stimulation of keratinocyte migration during cutaneous injury repair [61]. We also noticed that 6(OH)Mel can enhance the C-X-C motif chemokine ligand 1 (CXCL1) level ( Figure 4D). CXCL1 is an antimicrobial protein, which also plays a role in the antiviral innate immune response [62]. These findings are consistent with the above bioinformatic analyses and well-documented immunomodulatory and protective properties of melatonin [63].

Results of Telomerase Assay
The RNAseq data and pathway analysis (Table 4) were confirmed by measurement of telomerase activity by TRAP assay in HaCaT keratinocytes treated with melatonin and its metabolites for 24 h at a concentration of 10 −4 M ( Figure 5). The use of an immortalized keratinocyte line was dictated by its phenotypic stability during passaging. The influence of the melatonin on telomerase activity has also been indicated by GO analysis of RNAseq data after stimulation with melatonin, which included overrepresentation of Positive-Regulation-of-Telomerase-Activity gene set enrichment with the parameter of Rank at Max equal 4064. These data and bioinformatics analyses support previous studies showing antiaging activity of melatonin in the human skin [31,43].

Energy Yielding Metabolism Assays
The bioinformatics analyses presented in Sections 3.1 and 3.2 that indicated an effect on energy yielding metabolism were substantiated by qPCR analyses. Specifically, melatonin and its hydroxymetabolites stimulated sirtuin 1 (SIRT-1) and PGC1α, while inhibiting LDHA gene expressions in HaCaT keratinocytes ( Figure 6). The latter finding has been confirmed by inhibition of lactate production in cultured cells treated with 6(OH)Mel and AFMK ( Figure 6D).

Results of Telomerase Assay
The RNAseq data and pathway analysis (Table 4) were confirmed by measurement of telomerase activity by TRAP assay in HaCaT keratinocytes treated with melatonin and its metabolites for 24 h at a concentration of 10 −4 M ( Figure 5). The use of an immortalized keratinocyte line was dictated by its phenotypic stability during passaging. The influence of the melatonin on telomerase activity has also been indicated by GO analysis of RNAseq data after stimulation with melatonin, which included overrepresentation of Positive-Regulation-of-Telomerase-Activity gene set enrichment with the parameter of Rank at Max equal 4064. These data and bioinformatics analyses support previous studies showing antiaging activity of melatonin in the human skin [31,43].

Energy Yielding Metabolism Assays
The bioinformatics analyses presented in Sections 3.1 and 3.2 that indicated an effect on energy yielding metabolism were substantiated by qPCR analyses. Specifically, melatonin and its hydroxymetabolites stimulated sirtuin 1 (SIRT-1) and PGC1α, while inhibiting LDHA gene expressions in HaCaT keratinocytes ( Figure 6). The latter finding has been confirmed by inhibition of lactate production in cultured cells treated with 6(OH)Mel and AFMK ( Figure 6D).     Next, we performed a comprehensive metabolic analysis using the Seahorse metabolic analyzer platform. The obtained data on the oxygen consumption rate (mitochondrial respiration, OCR, Figure 6A) and extracellular acidification rate (glycolytic function, ECAR, Figure 6B) showed consistency with the RNAseq, bioinformatics and functional analyses presented above. We noticed decreased basal respiration and ATP production 24 h following AFMK or 6(OH)Mel compared to control or melatonin stimulated conditions ( Figure 7A). Maximal respiration, proton leak and spare respiratory capacity were unaffected by melatonin or its metabolites. Both AFMK and 6(OH)Mel, but not melatonin itself, significantly decreased glycolysis and glycolytic production in HEKn cells with a low oxygen consumption rate ( Figure 7B), which is consistent with inhibition of lactate production by the compounds presented in Figure 5E. In general, these studies are in line with recognized abilities of melatonin and its metabolites to affect mitochondrial functions and their intramitochondrial metabolism [2,9,[14][15][16]22,27,42,64].
Next, we performed a comprehensive metabolic analysis using the Seahorse metabolic analyzer platform. The obtained data on the oxygen consumption rate (mitochondrial respiration, OCR, Figure 6A) and extracellular acidification rate (glycolytic function, ECAR, Figure 6B) showed consistency with the RNAseq, bioinformatics and functional analyses presented above. We noticed decreased basal respiration and ATP production 24 h following AFMK or 6(OH)Mel compared to control or melatonin stimulated conditions ( Figure 7A). Maximal respiration, proton leak and spare respiratory capacity were unaffected by melatonin or its metabolites. Both AFMK and 6(OH)Mel, but not melatonin itself, significantly decreased glycolysis and glycolytic production in HEKn cells with a low oxygen consumption rate ( Figure 7B), which is consistent with inhibition of lactate production by the compounds presented in Figure 5E. In general, these studies are in line with recognized abilities of melatonin and its metabolites to affect mitochondrial functions and their intramitochondrial metabolism [2,9,[14][15][16]22,27,42,64]. While previous studies on metabolic functions were predominantly performed on isolated mitochondria, the Seahorse metabolic platform utilizes the intact cells. The reported inhibition of basal respiration and ATP production and of glycolysis by melatonin metabolites in keratinocytes is consistent with reported previously inhibition of cell proliferation by melatonin and its downstream metabolites [13,21,30,44]. However, the lack While previous studies on metabolic functions were predominantly performed on isolated mitochondria, the Seahorse metabolic platform utilizes the intact cells. The reported inhibition of basal respiration and ATP production and of glycolysis by melatonin metabolites in keratinocytes is consistent with reported previously inhibition of cell proliferation by melatonin and its downstream metabolites [13,21,30,44]. However, the lack of an effect of melatonin on the above metabolic parameters requires further investigation since melatonin is metabolized rapidly within mitochondria to 6(OH)Mel and AFMK and the presented cellular effect may represent an average of both melatonin effect and that of its metabolites. This issue will require additional studies that are beyond the scope of this manuscript.

Conclusions
Comprehensive analysis of RNAsequence data has shown significant changes in gene expression induced by melatonin and its two indolic (6(OH)Mel) and kynuric (AFMK) metabolites with an overlapping and differential pattern. This resulted in common regulation of top diseases and disorders and of molecular and cellular functions with some additional specificity for each of the tested molecules. This is common for melatonin, 6(OH)Mel and AFMF signaling pathways, with some signaling pathways being specific for each tested molecule. These may include common and separate nuclear receptors or regulatory proteins on which these molecules will act. GSEA, Reactome FIViZ and IPA indicated activation of p53 signaling and NRF2-mediated antioxidative pathways, which is consistent with our previous studies showing protective and antioxidative functions of melatonin and its metabolites. An overrepresentation of circadian clock and antiaging signaling by melatonin derivatives is consistent with antiaging properties of the melatonin. The GO showing an activation of the keratinocyte differentiation program by melatonin and the GSEA indicating antitumor activity are consistent with functional data reported previously demonstrating antiproliferative properties of melatonin and its metabolites in skin cells. Interestingly, the GSEA and IPA indicated, respectively, antilipidemic potential and antiviral responses controlled by melatonin and its metabolites. In addition, melatonin and metabolites had an effect on cellular metabolism and mitochondria functions in keratinocytes.
It is expected that the above changes will affect keratinocyte proliferation and differentiation programs to improve epidermal barrier formation, leading to increased protection against environmental insults and different pathogens. In addition, these analyses suggest possible applications of melatonin, AFMK and 6(OH)Mel in translational medicine.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/antiox10040618/s1, Table S1: Summary of Ingenuity Pathway Analysis (IPA) for top biological functions regulated by melatonin and its metabolites based on the data obtained from primary human epidermal keratinocytes. Figure S1: Volcano Plot diagram of differentially expressed genes for melatonin vs EtOH (A), 6-hydroxymelatonin vs EtOH (B), and AFMK vs EtOH (C).  Data Availability Statement: All raw RNAseq data were deposited in NCBI GEO (GSE147588).

Conflicts of Interest:
The authors declare no conflict of interest.