Abstract We analysed the whole-genome transcriptional profile of 6 cell lines of dark melanocytes (DM) and 6 of light melanocytes (LM) at basal conditions and after ultraviolet-B (UVB) radiation at different time points to investigate the mechanisms by which melanocytes protect human skin from the damaging effects of UVB. Further, we assessed the effect of different keratinocyte-conditioned media (KCM+ and KCM-) on melanocytes. Our results suggest that an interaction between ribosomal proteins and the P53 signaling pathway may occur in response to UVB in both DM and LM. We also observed that DM and LM show differentially expressed genes after irradiation, in particular at the first 6h after UVB. These are mainly associated with inflammatory reactions, cell survival or melanoma. Furthermore, the culture with KCM+ compared with KCM- had a noticeable effect on LM. This effect includes the activation of various signaling pathways such as the mTOR pathway, involved in the regulation of cell metabolism, growth, proliferation and survival. Finally, the comparison of the transcriptional profiles between LM and DM under basal conditions, and the application of natural selection tests in human populations allowed us to support the significant evolutionary role of MIF and ATP6V0B in the pigmentary phenotype. Introduction Melanocytes are melanin-producing cells that, in addition to hold a major role in the pigmentary phenotype, also play an important part in the protection of the skin against the damaging effects of ultraviolet-B (UVB) radiation, such as erythema, sunburn, development of malignant melanoma or other skin cancers [[41]1–[42]4]. The advent of cDNA microarray technology has allowed a preliminary understanding of the gene interactions and regulatory networks that take place in pigmentary cells in response to UVB [[43]5–[44]7]. One of the first reports using cDNA microarrays in various cell lines of human melanocytes for around 9,000 human genes [[45]5] showed that various genes, mainly related to DNA/RNA synthesis and modification, ribosomal proteins or solute carriers and ionic channels, were modulated 4 hours after a single dose of UVB irradiation (100mJ/cm^2). Later, Yang et al. [[46]6] using a higher density microarray (with probes for approximately 47,000 transcripts), although for a single cell line of melanocytes, analysed the response of melanocytes to UVB. In contrast to Valéry et al. [[47]5], Yang et al. [[48]6] selected a 24-hour time point after UVB irradiation and reported a set of p53-target genes as major agents involved in the UV response. However, many questions remain unsolved yet. For example, although the damage and the collateral consequences of UVB in the human skin are known to differ among individuals of different geographical origin and skin color [[49]8], the different transcriptional responses that could arise between cultured human melanocytes from dark and light donors (hereinafter DM and LM, respectively) have not been completely elucidated. A recent work [[50]9] performed a genome-wide transcriptome analysis of both DM and LM under basal conditions using RNA-Seq technology and found only 16 genes differentially expressed in the two cell types. However, their results could be somehow limited by the small number of melanocyte lines of each type (2 DM and 2 LM) analysed. Furthermore, the response of melanocytes to UV radiation is known to be mediated by paracrine factors released by keratinocytes, which modulate the growth rate and dendricity of melanocytes, and which ultimately lead to an increased production of melanin [[51]10–[52]18]. In some cases this has been shown by growing melanocytes with keratinocyte conditioned media (KCM) in vitro; however, the procedure by which this medium is obtained varies among studies. Thus, while in some experiments this medium is collected after the irradiation of keratinocytes (hereinafter KCM+) [[53]19–[54]20], in other studies the medium is collected from keratinocytes that have not been previously irradiated (hereinafter KCM-) [[55]10, [56]21–[57]22]. Given the current lack of a consensus to define how to collect this media, we aimed to analyse the putative different responses that might arise when culturing melanocytes with either KCM- or KCM+. Therefore, the objective of this work was to achieve a full view of the regulatory mechanisms that melanocytes undergo in response to UVB. Thus, we analyse herein the whole-genome transcriptional profile of dark and light melanocytes under basal conditions and after UVB irradiation at different time points (6, 12 and 24 hours) by means of gene expression microarrays. Further, we also aimed to assess the effect of different keratinocyte-conditioned media on melanocytes at a whole-genome level. With that aim, melanocytes were cultured in medium supplemented with keratinocyte-conditioned medium obtained both from non-irradiated (KCM-) and irradiated keratinocytes (KCM+). This work outperforms previous studies in many regards: 1) we interrogate a large number of probes in the genome, including genes (28,000) and other non-coding RNAs (7,419), 2) we include both DM and LM and assess their transcriptional differences, 3) importantly, we use a relatively high number of biological replicates (6 cell lines of DM and 6 of LM), which minimises the noise from variability among individuals, 4) we perform a time-series analysis that detects both early and later stress responses and 5) we cultivate melanocytes with KCM- and KCM+ and assess their distinct influence. Materials and Methods Cell cultures Human epidermal keratinocytes were purchased from Cascade Biologics (Life technologies, Carlsbad, CA, USA). Cells were cultured in EpiLife Medium supplemented with human keratinocyte growth supplement (HKGS). Human epidermal melanocytes were also purchased from Cascade Biologics: six lines isolated from lightly pigmented neonatal foreskin (LM), and six lines from darkly pigmented neonatal foreskin (DM). These melanocytes were cultivated in Medium 254 supplemented with 1% human melanocyte growth supplement (HMGS). All the cell lines were maintained in an incubator under an atmosphere of 5% CO[2] at 37°C. Media were refreshed every two days. UV irradiation and Keratinocyte-conditioned medium UV irradiation was performed in an ICH2 photoreactor (LuzChem, Canada) at 37°C. Cultures were irradiated at 75 mJ/cm^2 UVB, based on our previous work [[58]20], as we observed that this dosage led to a notable physiological effect but did not affect cell viability in both keratinocytes and melanocytes. Keratinocyte supernatants were harvested from both non-irradiated (KCM-) and irradiated keratinocytes (24 hours after treatment) (KCM+) and kept frozen at -80°C until subsequent use. Subconfluent melanocyte cultures were cultivated in Medium 254 supplemented with HMGS and KCM+ or KCM- medium in a proportion 1:1. The following day they were irradiated with 75 mJ/cm^2 of UVB, and harvested at 6, 12 and 24 hours post irradiation. We used non-irradiated control cultures that were covered by aluminium foil during irradiation ([59]Fig 1). Fig 1. Graphical scheme of the experimental design. [60]Fig 1 [61]Open in a new tab Microarrays RNA from irradiated and non-irradiated melanocytes was extracted using the RNA extraction kit from Ambion (Life technologies). Samples were quantified using a UV/VIS NanoDrop 8000 (Thermo Fisher, Waltham, MA, USA), and RNA integrity was analysed through an Agilent 2100 Bioanalyzer using Agilent RNA 6000 Nano Chips (Agilent Technologies, Santa Clara, CA, USA). For each labeling reaction 100 ng RNA were used, with the Low input Quick Amp Labeling kit, one color (Agilent Technologies). First, total RNA was retrotranscribed using AffinityScript Reverse Transcriptase (Agilent Technologies) and Oligo dT primers linked to promoter T7. The synthesized double stranded cDNA was in vitro transcribed by T7 RNA polymerase with Cy3-CTP in order to achieve labeled and amplified cRNA. These samples were purified with RNeasy Mini kit columns (Qiagen, Hilden, Germany) and quantified to determinate the yield (which should be higher than 0.825 μg per reaction) and the specific activity of the fluorochrome Cyanine 3 (which should be higher than 6 pmol/μg). All the samples satisfied these requirements. Samples were analysed using SurePrint G3 Human GE Microarrays (Agilent Technologies), which have probes for 27,958 annotated genes and 7,419 long intergenic non-coding RNAs (lincRNAs). The hybridization step was performed using the SureHyb hybridization chamber (Agilent Technologies) and 600 ng of labeled cRNA samples, for 17 hours at 65°C and 10,000 rpms in a hybridization oven. Microarrays were stabilized with ozone-barrier slide covers (Agilent Technologies). Image processing of the microarrays was performed by using the Agilent Feature Extraction software v10.7.3.1. This software performs 9 evaluation parameters to check the quality of the microarrays. The quality control parameters included, among others, the coefficient of variation of the processed signal from non-control probes and spike-ins (%CV), the percentage of outlier probes as regards the replicated probes population, the intensity of the signal of the negative controls and the limit of detection and linearity of the Spike-Ins signal. Microarray data pre-processing and normalization Raw data were processed with GeneSpring GX software v11.5.1 (Agilent Technologies). Feature extraction flags were transformed as follows: if feature was not positive and significant, not uniform, not well above background or was a population outlier: compromised; if feature was saturated: not detected. We performed a variance-stabilizing transformation of the data, which is a key step, but often not considered, in the pre-processing of microarrays data. Most of the subsequent statistical analyses assume that the data follow a normal distribution, with a constant variance independent of the mean of the data. Gene-expression microarray data, however, often have a variance that changes non-linearly with the mean, and thus, log transformations, which are used in the transformation of these data, can inflate the variance of observations near the background. Thus, our data were subjected to a DDHF (Data-Driven Haar-Fisz) transformation for variance stabilization with the R package DDHFm [[62]23]. This method stabilizes the variance of replicated intensities from microarray data and produces transformed intensities that are much closer to the Gaussian distribution than other methods. Furthermore, it can be adapted to different or uncertain distributions, and therefore, it is ideal for the variance stabilization of microarray data. Data were transformed to log base 2 and normalized following the quantile method [[63]24]. Flag spot information in data files was used to filter probe sets. Entities in which more than 50% of samples in 1 out of any 7 conditions (0h, 6h KCM-, 12h KCM-, 24h KCM, 6h KCM+, 12h KCM+ and 24h KCM+) had “detected” flags were maintained for the analyses. Quality (QC) Metrics and Principal Component Analysis (PCA) QC-Metrics was performed with GeneSpring GX software. Gene expression of the transformed and normalized data were subjected to unsupervised classification by means of Principal Component Analysis (PCA) as a preliminary exploratory approach to detect outliers, or the existence of defined clusters based on time points, pigmentation of the cells or the type of KCM used for culture. We used The Unscrambler X v10.3 (CAMO A/S, Trondheim, Norway) and applied the full cross validation method to estimate the stability and performance of the model. Comparison of expression profiles Statistical analysis for the comparison of expression profiles was performed with SAM (Significance Analysis of Microarrays, [[64]25]), using two class non-pairwise comparisons and 500 permutations in each test. The significance of the tests was given by the lowest False Discovery Rate at which the gene is called significant based on [[65]26], adjusted for multiple tests. Pathway enrichment analysis Enrichment analysis was performed using Web-based Gene Set Analysis Toolkit (WebGestalt) ([66]http://bioinfo.vanderbilt.edu/-webgestalt/option.php), using all the probes analyzed in the microarray as the reference list, and The Kyoto Encyclopedia of Genes and Genomes (KEGG) database of pathways. The significance analysis was performed using the Hypergeometric test. P-values were corrected for multiple tests following the Bonferroni procedure. The minimum number of genes for enrichment was set at 5, and the significance level at Bonferroni adjusted-p<0.01, in order to be conservative, avoid false positives and achieve more robust results. Validation by RT-qPCR We selected 6 genes showing a change in expression between dark and light melanocytes or after UV irradiation in the microarrays for validation with Real-Time quantitative PCR (RT-qPCR). cDNA was synthesized from 2 μg of total extracted RNA using the First Strand cDNA Synthesis Kit (ThermoFisher) and was used as a template for RT-qPCR analyses. Four different cell lines were analysed (2 of dark melanocytes, and 2 of light melanocytes). RT-qPCR reactions were performed with SYBR Green in a StepOne thermocycler (Life Technologies). Primer sequences (5'-3') were the following: MIFf_GAAGGCCATGAGCTGGTCT, MIFr_GGTTCCTCTCCGAGCTCAC, FDXRf_CTGAGGCAGAGTCGAGTGAAG, FDXRr_CCCGAAGCTCCTTAATGGTGA, TP53I3f_AATGCTTTCACGGAGCAAATTC, TP53I3r_TTCGGTCACTGGGTAGATTCT, ATP6VOBf_CCATCGGAACTACCATGCAGG, ATP6VOBr_TCCACAGAAGAGGTTAGACAGG, MDM2f_GAATCATCGGACTCAGGTACATC, MDM2r_TCTGTCTCACTAATTGCTCTCCT, RPL6f_ATTCCCGATCTGCCATGTATTC and RPL6r_TACCGCCGTTCTTGTCACC. Thermocycling conditions were optimized for each pair of primers to obtain 95–100% efficiency and r^2>0.99 in the reaction. Gene expression was normalized to the housekeeping gene GAPDH. Each reaction was performed in triplicate and values were averaged to calculate relative expression levels. Selection tests Preliminary screenings to detect deviations from neutrality were performed using the 1000 Genomes Selection Browser ([67]http://hsb.upf.edu/) [[68]27], which implements several neutrality tests (Tajima’s D, Fay & Wu’s H, Fu’s F, Fu and Li’s F*, Fu and Li’s D* and EHH, among many others) and provides genome based rank “p-values”, that help to identify which SNPs or regions have significantly high scores compared to the rest of the genome. Further selection tests in candidate loci were performed with DnaSP [[69]28]. We obtained the genotypes of the European (n = 760 chromosomes), African (n = 492) and Asian (n = 572) populations from the 1000 Genomes Project (Phase I May 2011) using SPSmart v5.1.1 (dbSNP build 132) [[70]29]. The orthologous sequence of the chimpanzee was obtained from the UCSC Genome Browser and aligned to the human sequences with ClustalW. For each population, we calculated Tajima’s D, Fu & Li’s D and Fay & Wu’s H with DnaSP [[71]28]. P-values for these tests were obtained using the interface for standard coalescent simulations conditioned on the number of segregating sites. Results and Discussion Quality metrics and PCA QC-Metrics revealed 2 outlier arrays that did not satisfy the quality parameters: L_5.6K- (LM; replicate_5; 6h; KCM-) and L_4.24K- (LM; replicate_4; 24h; KCM-). Thus, those samples were removed from the subsequent statistical analyses. Second, we performed a Principal Component Analysis (PCA), an exploratory multivariate statistical technique for simplifying complex data sets [[72]30], that has been used for the analysis of microarray data in search of outlier genes [[73]31] or to identify temporal patterns in time-series analyses [[74]32]. The PCA ([75]Fig 2) allowed us to have an overview of the temporal patterns or differentially expressed genes between dark and light melanocytes, or between the culture with KCM+ or KCM-. It showed an apparent general homogeneity, revealing no additional potential outliers and a coherent clustering of our samples according to different variables, which was valuable to discard the presence of outliers or experimental errors. The PCA showed a time-point clustering defined by the second component, revealing a major separation of the samples at 6 hours, while the samples corresponding to 12 and 24 hours clustered close to the controls (0 hours), thus suggesting an early response from melanocytes to UVB that returned again to basal levels after the first 6 hours. At 6 hours we also observed a differential response according to pigmentation defined by the first component. At 24 hours, an apparent clustering regarding the culture of light melanocytes with KCM- or KCM+ was also noticed. Fig 2. Principal Component Analysis. Fig 2 [76]Open in a new tab Charts a) and b) show the same 2-dimesional representation of the data according to the first 2 principal components, but colored according to different variables. Thus, in a) the effects of time (Squares: Time = 0; Dots: Time = 6 hours; Triangles: Time = 12 hours; Diamonds: Time = 24 hours) and pigmentation (Yellow = Light melanocytes; Brown = Dark melanocytes) are highlighted, while in b), it is the time (Squares: Time = 0; Dots: Time = 6 hours; Triangles: Time = 12 hours; Diamonds: Time = 24 hours) and the type of KCM used which are highlighted (Green: KCM-; Red: KCM+). Identification of differentially expressed probes after UVB A total of 26,493 probes were examined per microarray. By means of SAM [[77]25] we identified the statistically significant differentially expressed genes. Because probes may correspond to both genes and non-coding RNAs, we explicitly indicated when they corresponded to non-coding RNA. We first looked for common genes differentially expressed in DM and LM across time after UVB irradiation. We focused on the top upregulated and downregulated genes at 6, 12 and 24h, and in order to provide robust results, we identified those genes that were significantly up- or downregulated at more than one point (Tables [78]1 and [79]2). The adjusted p-value for all these genes was <0.0001. Table 1. Most significantly upregulated genes at more than one time point (non-coding RNAs are indicated with *) (Bonferroni-adjusted p-value <0.0001). Time points Gene symbol Accession number Description 6, 12 and 24h FDXR [80]NM_004110 ferredoxin reductase, nuclear gene encoding mitochondrial protein EPHA2 [81]NM_004431 EPH receptor A2 RPL6 [82]NM_001024662 ribosomal protein L6 VWCE [83]NM_152718 von Willebrand factor C and EGF domains UBD [84]NM_006398 ubiquitin D CXCL1 [85]NM_001511 chemokine (C-X-C motif) ligand 1 (melanoma growth stimulating activity, alpha) MDM2 [86]NM_002392 Mdm2 p53 binding protein homolog (mouse), transcript variant MDM2 RPL41 [87]NM_001035267 ribosomal protein L41 TNFRSF10C [88]NM_003841 tumor necrosis factor receptor superfamily, member 10c DDB2 [89]NM_000107 damage-specific DNA binding protein 2, 48kDa GRM2 [90]NM_000839 glutamate receptor, metabotropic 2 TP53I3 [91]NM_004881 tumor protein p53 inducible protein 3 ISCU [92]NM_014301 iron-sulfur cluster scaffold homolog (E. coli) 6 and 12h GADD45A [93]NM_001924 growth arrest and DNA-damage-inducible, alpha PLK3 [94]NM_004073 polo-like kinase 3 BTG2 [95]NM_006763 BTG family, member 2 TRIAP1 [96]NM_016399 TP53 regulated inhibitor of apoptosis 1 PIDD [97]NM_145886 p53-induced death domain protein CDKN1A [98]NM_078467 cyclin-dependent kinase inhibitor 1A (p21, Cip1) TP53INP1 [99]NM_033285 tumor protein p53 inducible nuclear protein 1 SESN1 THC2527965 Sestrin 1, partial (68%) BAG1 [100]NM_004323 BCL2-associated athanogene XPC [101]NM_004628 xeroderma pigmentosum, complementation group C *lincRNA chrX:64042150–64093950 12 and 24h RPS2 [102]NM_002952 ribosomal protein S2 RPL26 THC2550570 ribosomal protein L26, partial (91%) PLXNB2 [103]NM_012401 plexin B2 KRT17 [104]NM_000422 keratin 17 ACTA2 [105]NM_001613 actin, alpha 2, smooth muscle, aorta SULF2 [106]NM_018837 sulfatase 2 PVRL4 [107]NM_030916 poliovirus receptor-related 4 CSRP2 [108]NM_001321 cysteine and glycine-rich protein 2 DRAM1 [109]NM_018370 DNA-damage regulated autophagy modulator 1 BBC3 [110]NM_014417 BCL2 binding component 3 *LOC344887 [111]NR_033752 NmrA-like family domain containing 1 pseudogene, non-coding RNA RELB [112]NM_006509 v-rel reticuloendotheliosis viral oncogene homolog B *LOC642335 [113]AK098072 cDNA FLJ40753 fis, clone TRACH2001188. KIAA1324 [114]NM_020775 KIAA1324 NOV [115]NM_002514 nephroblastoma overexpressed gene RPS27L [116]NM_015920 ribosomal protein S27-like PRODH [117]NM_016335 proline dehydrogenase (oxidase) 1, nuclear gene encoding mitochondrial protein 6 and 24h GDF15 [118]NM_004864 growth differentiation factor 15 NFKBIA [119]NM_020529 nuclear factor of kappa light polypeptide gene enhancer in B-cells inhibitor, alpha [120]Open in a new tab Table 2. Most significantly downregulated genes at more than one time point (non-coding RNAs are indicated with *) (Bonferroni-adjusted p-value <0.0001). Time points Gene symbol Accession number Description 6, 12 and 24h LGALS3 [121]NM_002306 lectin, galactoside-binding, soluble, 3 PDSS2 [122]NM_020381 prenyl (decaprenyl) diphosphate synthase, subunit 2 MAGI3 [123]NM_152900 membrane associated guanylate kinase, WW and PDZ domain containing 3 PSD3 [124]NM_015310 pleckstrin and Sec7 domain containing 3 *lincRNA chr10:114583921–114587485 forward strand 6 and 12h SBF2 [125]NM_030962 SET binding factor 2 XRCC4 [126]NM_022550 X-ray repair complementing defective repair in Chinese hamster cells 4 VAV2 [127]NM_003371 vav 2 guanine nucleotide exchange factor ROR1 [128]NM_005012 receptor tyrosine kinase-like orphan receptor 1 ERC1 [129]NM_178040 ELKS/RAB6-interacting/CAST family member 1 *lincRNA chr17:67547498–67549996 forward strand 12 and 24h HMG20B [130]NM_006339 high mobility group 20B TUBA1B [131]NM_006082 tubulin, alpha 1b SMYD3 [132]NM_022743 SET and MYND domain containing 3 SCFD2 [133]NM_152540 sec1 family domain containing 2 VTI1A [134]NM_145206 vesicle transport through interaction with t-SNAREs homolog 1A (yeast) PRR4 [135]NM_007244 proline rich 4 (lacrimal) PARD3 [136]NM_019619 par-3 partitioning defective 3 homolog (C. elegans) RABGAP1L [137]NM_014857 RAB GTPase activating protein 1-like TTC28 [138]NM_001145418 tetratricopeptide repeat domain 28 PCCA [139]NM_000282 propionyl CoA carboxylase, alpha polypeptide MAN1C1 [140]NM_020379 mannosidase, alpha, class 1C, member 1 A4GALT [141]NM_017436 alpha 1,4-galactosyltransferase MSRA [142]NM_012331 methionine sulfoxide reductase A ANO4 [143]NM_178826 anoctamin 4 SSBP2 [144]NM_012446 single-stranded DNA binding protein 2 STX8 [145]NM_004853 syntaxin 8 REXO1 [146]NM_020695 REX1, RNA exonuclease 1 homolog (S. cerevisiae) SH3KBP1 [147]NM_031892 SH3-domain kinase binding protein 1 BBS9 [148]NM_198428 Bardet-Biedl syndrome 9 BCKDHB [149]NM_000056 branched chain keto acid dehydrogenase E1, beta polypeptide SORCS1 [150]NM_001206572 sortilin-related VPS10 domain containing receptor 1 TPK1 [151]NM_022445 thiamin pyrophosphokinase 1 LINGO2 [152]NM_152570 leucine rich repeat and Ig domain containing 2 FRY [153]NM_023037 furry homolog (Drosophila) PDE3B [154]NM_000922 phosphodiesterase 3B, cGMP-inhibited KCNQ2 [155]NM_172109 potassium voltage-gated channel, KQT-like subfamily, member 2 PPIA THC2525667 Peptidylprolyl isomerase A *lincRNA chr2:7214634–7218011 reverse strand *lincRNA chr4:79892901–80229698 forward strand *lincRNA chr18:42263052–42278652 forward strand *lincRNA chr18:74178337–74203637 forward strand *lincRNA chr7:125564239–125734564 forward strand 6 and 24h FAM78B [156]NM_001017961 family with sequence similarity 78, member B VAV3 [157]NM_006113 vav 3 guanine nucleotide exchange factor [158]Open in a new tab Common upregulated genes after UVB irradiation Some of the genes included in this category ([159]Table 1) have already been reported to be associated with the response to ultraviolet irradiation, which gives robustness to our inferences. The most significantly upregulated gene was FDXR, which serves as the first electron transfer protein in all the mitochondrial P450 systems, and has been reported to be upregulated in response to UV irradiation damage in dendritic cells [[160]33] and melanocytes [[161]6]. Importantly, we also observed several genes involved in the regulation of the cell cycle, in the response to stress, in the repair of DNA damage caused by UV that can lead to xeroderma pigmentosum, as well as genes that are associated with melanoma. We also observed several genes that take part in the regulation of the cell cycle and in the cellular response to stress and that are directly or indirectly involved in the p53 pathway. Some of them modulate P53-mediated apoptosis or cell death in response to stresses like UV irradiation or DNA damage, like TP53I3, PLK3, TRIAP1, PIDD, CDKN1A, TP53INP1, SESN1, BBC3, TNFRSF10C, DRAM1 and MDM2. Other genes that also are upregulated and participate in UV irradiation-induced apoptosis include RELB and EPHA2. Another group of the upregulated genes are components of the nucleotide excision repair (NER) pathway that are associated with the reparation of DNA damage caused by UV, and which include XPC or DDB2. Malfunction of these genes can lead to xeroderma pigmentosum, a recessive disease that is characterized by an increased sensitivity to UV light and a high predisposition for skin cancer development. Several other genes among the top upregulated ones have been reported to be directly or indirectly associated with melanoma, such as BTG2, BAG1, CXCL1, PLXNB2, CSRP2, PRODH or GDF15. Various genes that encode ribosomal proteins such as RPL6, RPL41, RPS2, RPL26 and RPS27L were also observed. Intriguingly, we observed an upregulation of the gene NOV. The protein encoded by this gene is of particular relevance as it has been reported to be essential for the correct development and growth of melanocytes [[162]34]. During development, melanocytes migrate to the epidermis and attach to the basement membrane upon contact with keratinocytes. Development of melanocytes must be tightly regulated and must remain at stable levels in relation to keratinocytes. Fukunaga-Kalabis et al. [[163]34] discovered that NOV is upregulated in melanocytes upon contact with keratinocytes in culture, mediating the growth inhibition of melanocytes in order to regulate their spatial location and number. Our results suggest that this gene could also participate in the regulation of melanocytes' growth in response to UVB, most likely by inhibiting their proliferation and allowing either the triggering of cell death or reparation events. Common downregulated genes after UVB irradiation Among the top downregulated genes in response to UVB irradiation ([164]Table 2), of particular interest was LGALS3, which plays a role in numerous cellular functions including apoptosis, innate immunity, cell adhesion and T-cell regulation, and regulates the expression of several genes that are aberrantly expressed in highly aggressive melanoma cells [[165]35]. Another interesting downregulated gene was PDSS2, which encodes the prenyl side-chain of coenzyme Q (CoQ), one of the key elements in the respiratory chain. As it has been reported that UV light depletes CoQ10 from the skin [[166]36], this consequently suggests that the downregulation of PDSS2 could be one of the first inducers of reactive oxygen species (ROS) production and, consequently, of oxidative damage to the DNA in the cells ultimately caused by UVB. Several neuron-related genes were also downregulated by UVB, like ROR1, ERC1, PARD3, SORCS1, LINGO2 and KCNQ2. As neurons and melanocytes share the same embryonic origin (the neural crest) they might likely share some cell regulatory processes [[167]37]. Our results suggest that some genes that are involved in neuronal growth or migration could also be found in melanocytes exerting similar functions. In this regard, it was noticeable that a handful of lincRNAs were also downregulated after UVB irradiation. Although for most lincRNAs biological functions and mechanisms of action remain unknown, our results suggest that some lincRNAs are likely key elements of the regulatory machinery of melanocytes. Differential transcriptional profile of dark and light melanocytes 6 hours after UVB As inferred from the unsupervised PCA ([168]Fig 2) the greatest differences among melanocytes were found at 6 hours after UVB irradiation. From that point, melanocytes seem to have started to go back to the basal state. Thus, we focused on determining the differential expression between DM and LM at 6 hours after irradiation (the results for other time points can be found in [169]S1–[170]S4 Tables). Upregulated genes in LM vs DM 6 hours after UVB irradiation The significant upregulation of EDA2R in LM ([171]Table 3) suggests a putative role for this gene in response to UVB irradiation in light skinned individuals. EDA2R has been reported to be affected by recent natural positive selection [[172]38–[173]39], and its paralog, EDAR, has been strongly associated with skin pigmentation variability in humans [[174]40]. Strikingly, the intergenic region between EDA2R and the next gene on the same chromosome (AR) is the most divergent genomic segment between Africans and East Asians in the human genome [[175]41]. Table 3. Top 50 upregulated genes in LM vs DM at 6 hours after UVB irradiation (non-coding RNAs are indicated with *) (Bonferroni-adjusted p-value <0.0001). Gene symbol Accession number Description CDKN2A [176]NM_000077 cyclin-dependent kinase inhibitor 2A (melanoma, p16, inhibits CDK4) SNAR-A3 [177]NR_024214 small ILF3/NF90-associated RNA A3, small nuclear RNA KIAA1377 [178]NM_020802 uncharacterized protein KIAA1377 TTC18 [179]NM_145170 tetratricopeptide repeat domain 18 NCAM1 [180]NM_001242607 neural cell adhesion molecule 1, tr. variant 5 HOXB13 [181]NM_006361 homeobox B13 PYCARD [182]NM_013258 PYD and CARD domain containing CSRNP3 [183]NM_001172173 cysteine-serine-rich nuclear protein 3 PKMYT1 [184]NM_182687 protein kinase, membrane associated tyrosine/threonine 1 *lincRNA:chr1:85932812–85974062 reverse strand QPCT [185]NM_012413 glutaminyl-peptide cyclotransferase EDA2R [186]NM_001242310 ectodysplasin A2 receptor SGMS1 [187]NM_147156 sphingomyelin synthase 1 RPL37A [188]NM_000998 ribosomal protein L37a HIST1H4L [189]NM_003546 histone cluster 1, H4l GDPD1 [190]NM_182569 glycerophosphodiester phosphodiesterase domain containing 1 HIST1H3B [191]NM_003537 histone cluster 1, H3b *lincRNA:chr10:133738235–133744210 forward strand GDPD5 [192]NM_030792 Glycerophosphodiester phosphodiesterase domain containing 5 SUV420H2 [193]NM_032701 suppressor of variegation 4–20 homolog 2 (Drosophila) MOBKL2B [194]NM_024761 MOB1, Mps One Binder kinase activator-like 2B (yeast) MXD3 [195]NM_031300 MAX dimerization protein 3 TUBA1B [196]NM_006082 tubulin, alpha 1b SAC3D1 [197]NM_013299 SAC3 domain containing 1 *LOC390595 [198]NM_001163692 ubiquitin-associated protein 1-like SNORD15A [199]NR_000005 small nucleolar RNA, C/D box 15A, small nucleolar RNA ZNF711 [200]NM_021998 zinc finger protein 711 *lincRNA:chrX:102139220–102156619 forward strand LTBP3 ENST00000525443 latent transforming growth factor beta binding protein 3 FAM164A [201]NM_016010 family with sequence similarity 164, member A ARHGEF10 ENST00000523711 Rho guanine nucleotide exchange factor (GEF) 10 S100B [202]NM_006272 S100 calcium binding protein B HMGN2 [203]NM_005517 high mobility group nucleosomal binding domain 2 C1orf15-NBL1 [204]NM_001204088 C1ORF15-NBL1 readthrough GALNT14 [205]NM_024572 polypeptide N-acetylgalactosaminyltransferase 14 (GalNAc-T14) SPTLC3 [206]NM_018327 serine palmitoyltransferase, long chain base subunit 3 IFI27L2 [207]NM_032036 interferon, alpha-inducible protein 27-like 2 RNF6 [208]NM_005977 ring finger protein (C3H2C3 type) 6 TUBB8 [209]NM_177987 tubulin, beta 8 PDGFRL [210]NM_006207 platelet-derived growth factor receptor-like ARPC5 ENST00000367534 actin related protein 2/3 complex, subunit 5, 16kDa PTGDS [211]NM_000954 prostaglandin D2 synthase 21kDa (brain) SLC2A13 [212]NM_052885 solute carrier family 2 (facilitated glucose transporter), member 13 CTSF [213]NM_003793 Cathepsin F *C1orf133 [214]NR_024337 SERTAD4 antisense RNA 1 WFDC1 [215]NM_021197 WAP four-disulfide core domain 1 TUBG1 [216]NM_001070 tubulin, gamma 1 SLC12A8 [217]NM_024628 solute carrier family 12, member 8 CXADR [218]NM_001338 coxsackie virus and adenovirus receptor SOX5 [219]NM_152989 SRY (sex determining region Y)-box 5 [220]Open in a new tab We also identified as upregulated some genes related to melanoma, like CDKN2A, whose expression has already been reported to be induced by UV radiation [[221]42] and which could be conferring light skinned individuals a higher susceptibility to develop melanoma in response to UV radiation, as well as several neuron-related genes and genes associated with the formation of tubulin, the major constituent of microtubules of the cytoskeleton, and which have been shown to mediate the transport the melanosomes inside the cell [[222]43]. Upregulated genes in DM vs LM 6 hours after UVB irradiation Next, we focused on the most significant upregulated genes in DM vs LM ([223]Table 4). In this case, we found several genes involved in inflammatory reactions. Some of them have been reported to be particularly involved in sunburn or inflammatory skin reactions in response to UVB, like IL6 [[224]44], PTGS2 [[225]45] or CCL2 [[226]46]. Similarly to LM, DM also showed an upregulation of various genes involved in melanoma progression as well as several genes related to the development of the central nervous system and neuronal processes. Table 4. Top 50 upregulated genes in DM vs LM at 6 hours after UVB irradiation (non-coding RNAs are indicated with *) (Bonferroni-adjusted p-value <0.0001). Gene symbol Accession numer Description HUMRPL26X THC2550570 ribosomal protein L26 partial (91%) MMP1 [227]NM_002421 matrix metallopeptidase 1 (interstitial collagenase) RPL7A [228]NM_000972 ribosomal protein L7a CCL2 [229]NM_002982 chemokine (C-C motif) ligand 2 NPTX2 [230]NM_002523 neuronal pentraxin II COL6A2 [231]NM_058174 collagen, type VI, alpha 2 S100A4 [232]NM_002961 S100 calcium binding protein A4 CXCL1 [233]NM_001511 chemokine (C-X-C motif) ligand 1 (melanoma growth stimulating activity, alpha) IL6 [234]NM_000600 interleukin 6 (interferon, beta 2) TMEM158 [235]NM_015444 transmembrane protein 158 (gene/pseudogene) PAMR1 [236]NM_015430 peptidase domain containing associated with muscle regeneration 1 TMEM8B [237]NM_001042590 collagen, type VI, alpha 2 CXCL5 [238]NM_002994 chemokine (C-X-C motif) ligand 5 TDRD9 [239]NM_153046 tudor domain containing 9 ANGPTL4 [240]NM_139314 angiopoietin-like 4 PMEPA1 [241]NM_020182 prostate transmembrane protein, androgen induced 1 *MEG3 [242]NR_003531 maternally expressed 3 (non-protein coding) non-coding RNA *C1D ENST00000412019 C1D nuclear receptor corepressor pseudogene C1QBP [243]NM_001212 complement component 1, q subcomponent binding protein *FAM27A [244]NR_024060 family with sequence similarity 27, member A, non-coding RNA TMEM132A [245]NM_017870 transmembrane protein 132A SLC16A2 [246]NM_006517 solute carrier family 16, member 2 (monocarboxylic acid transporter 8) ANPEP [247]NM_001150 alanyl (membrane) aminopeptidase *lincRNA:chr2:239460050–239536125 forward strand FN1 [248]NM_054034 fibronectin 1 MYOF [249]NM_133337 myoferlin NR4A3 [250]NM_173200 nuclear receptor subfamily 4, group A, member 3 EFS [251]NM_005864 embryonal Fyn-associated substrate GRTP1 [252]NM_024719 growth hormone regulated TBC protein 1 TNFRSF11B [253]NM_002546 tumor necrosis factor receptor superfamily, member 11b DEF6 [254]NM_022047 differentially expressed in FDCP 6 homolog (mouse) *FAM27A [255]NR_024060 family with sequence similarity 27, member A, non-coding RNA PTGS2 [256]NM_000963 prostaglandin-endoperoxide synthase 2 GUCA1B [257]NM_002098 guanylate cyclase activator 1B (retina) IL1RAP [258]NM_134470 interleukin 1 receptor accessory protein SUSD3 [259]NM_145006 sushi domain containing 3 *lincRNA:chr17:73585552–73590170 forward strand TNNI3 [260]NM_000363 troponin I type 3 (cardiac) C10orf116 [261]NM_006829 chromosome 10 open reading frame 116 SPON2 [262]NM_012445 spondin 2, extracellular matrix protein LYPD1 [263]NM_144586 LY6/PLAUR domain containing 1 IL27RA [264]NM_004843 interleukin 27 receptor, alpha FUT1 [265]NM_000148 fucosyltransferase 1 (galactoside 2-alpha-L-fucosyltransferase, H blood group) CARD9 [266]NM_052813 caspase recruitment domain family, member 9 1 SLC22A17 [267]NM_016609 solute carrier family 22, member 17 IL11 [268]NM_000641 interleukin 11 TNFAIP2 [269]NM_006291 tumor necrosis factor, alpha-induced protein 2 C15orf48 [270]NM_032413 chromosome 15 open reading frame 48 NT5E [271]NM_002526 5'-nucleotidase, ecto (CD73) CXCL3 [272]NM_002090 chemokine (C-X-C motif) ligand 3 [273]Open in a new tab An interesting observation was the upregulation of the lincRNA MEG3. The expression of this lincRNA, stimulated by cyclic-AMP (cAMP), seems to act as a growth suppressor in tumour cells through the activation of P53 [[274]47]. As UVR is one of the main stimulatory sources of cAMP, these results suggest that in response to UV radiation, DM upregulate the expression of MEG3 via cAMP liberation, which could confer protection against melanoma. Pathway enrichment analysis Focusing on single loci allows deciphering the differentially expressed genes between different categories (i.e. time or pigmentation). However, although this is useful to determine which genes can be key in the response to UVB, the full biological mechanisms underlying this response may remain obscure. Therefore, we used WebGestalt to look for pathways in KEGG (Kyoto Encyclopedia of Genes and Genomes) that were differentially overrepresented at each time point ([275]Table 5) in LM and DM. We observed that the most significant pathways overrepresented among the upregulated genes corresponded to ribosome and P53 signaling pathway in both LM and DM. Further analyses using other databases of pathways implemented in WebGestalt (Pathway Commons and Wikipathways) confirmed the involvement of these two pathways in the response to UVB (data not shown). Table 5. KEGG pathway enrichment analysis of the upregulated and downregulated genes in DM and LM after UVB irradiation. DM LM Pathway Adj-p Pathway Adj-p Upregulated at 6h p53 signaling pathway 9.49E-07 p53 signaling pathway 7.04E-07 Ribosome 3.00E-04 Ribosome 2.43E-06 Upregulated at 12h Ribosome 3.22E-11 Systemic lupus erythematosus 1.16E-06 RNA transport 7.11E-07 Ribosome 1.47E-06 Ribosome biogenesis in eukaryotes 7.00E-04 p53 signaling pathway 1.10E-03 p53 signaling pathway 2.13E-02 Mismatch repair 5.40E-03 Pathways in cancer 6.90E-03 Apoptosis 7.90E-03 Upregulated at 24h Ribosome 1.06E-09 Ribosome 7.38E-06 p53 signaling pathway 1.50E-03 Downregulated at 6h Adherens junction 1.51E-08 Ubiquitin mediated proteolysis 6.84E-12 Ubiquitin mediated proteolysis 7.21E-08 Adherens junction 5.23E-11 Wnt signaling pathway 9.12E-06 Endocytosis 6.14E-06 Colorectal cancer 3.53E-05 Wnt signaling pathway 3.14E-05 Progesterone-mediated oocyte maturation 4.04E-05 Progesterone-mediated oocyte maturation 5.32E-05 Endocytosis 6.24E-05 Cell cycle 7.36E-05 Pathways in cancer 7.60E-05 Insulin signaling pathway 5.00E-04 Systemic lupus erythematosus 1.00E-04 Oocyte meiosis 7.00E-04 Cell cycle 9.00E-04 Pathways in cancer 7.00E-04 Oocyte meiosis 1.00E-03 Colorectal cancer 1.20E-03 Endometrial cancer 1.10E-03 Neurotrophin signaling pathway 2.60E-03 Fc epsilon RI signaling pathway 1.70E-03 Fc gamma R-mediated phagocytosis 8.60E-03 B cell receptor signaling pathway 1.80E-03 Endometrial cancer 1.18E-02 ErbB signaling pathway 1.90E-03 Chronic myeloid leukemia 1.54E-02 Fc gamma R-mediated phagocytosis 2.20E-03 Bacterial invasion of epithelial cells 1.76E-02 T cell receptor signaling pathway 2.50E-03 Fc epsilon RI signaling pathway 2.11E-02 Insulin signaling pathway 5.30E-03 Chemokine signaling pathway 4.22E-02 Neurotrophin signaling pathway 1.97E-02 ErbB signaling pathway 4.22E-02 T cell receptor signaling pathway 4.22E-02 Downregulated at 12h Metabolic pathways 4.43E-05 Protein processing in endoplasmic reticulum 2.00E-03 Adherens junction 1.00E-04 Fc gamma R-mediated phagocytosis 3.00E-04 Propanoate metabolism 1.10E-03 Protein processing in endoplasmic reticulum 2.70E-03 Systemic lupus erythematosus 3.10E-03 Purine metabolism 4.90E-03 Regulation of actin cytoskeleton 1.20E-02 Valine, leucine and isoleucine degradation 4.30E-02 Pyrimidine metabolism 4.30E-02 Downregulated at 24h Cell cycle 2.20E-08 Systemic lupus erythematosus 8.56E-06 - DNA replication 2.54E-05 Oocyte meiosis 2.00E-03 Lysine degradation 2.12E-02 [276]Open in a new tab The role of P53, a tumour suppressor that promotes either cell cycle arrest and DNA repair, apoptosis or senescence [[277]48] in the response to UVB has already been reported [[278]6]. Our results are consistent with the proposed mechanism of P53 pathway regulation by ribosomal proteins [[279]49–[280]51]. Thus, we propose that under stress, there is an upregulation of the ribosomal biogenesis leading to an excess of ribosomal proteins that do not participate in the assembly of ribosomes. Instead, these translocate to the nucleoplasm where they interact with MDM2. Under normal conditions, MDM2 binds to the tumour suppressor P53 inhibiting its transcription. But if ribosomal proteins bind to MDM2, then the inhibition of P53 exerted by MDM2 is suppressed. The upregulation of MDM2 is usually modulated by P53 after the activation of P53-dependent targets, in order to inhibit the activity of P53 and thus restore the normal growth of the cell. However, if the stressing conditions are not completely restablished or DNA damage still exists in the cell, ribosomal proteins could continue interacting with MDM2 to allow to maintain the expression of P53 ([281]Fig 3). Among the ribosomal proteins that can bind to MDM2 are RPL5 [[282]52] and RPL11 [[283]53], both of which were among the upregulated ribosomal genes in this work. Fig 3. Proposed mechanism for the involvement of ribosomal proteins, MDM2 and p53 signaling pathway in the response to UVB irradiation. Fig 3 [284]Open in a new tab On the other hand, several pathways were significantly overrepresented among the genes that were downregulated after UVB exposure, especially in the first 6h in both DM and LM ([285]Table 5). Interestingly, the adherens junction pathway was downregulated in both cell types and at different time points. Adherens junctions play an important role maintaining skin homeostasis by mediating the interaction of melanocytes and keratinocytes, which control the proliferation of melanocytes [[286]54], thus preventing the development and progression of melanoma [[287]55]. The effect of keratinocyte conditioned medium Next, we assessed the expression profiles of melanocytes supplemented with keratinocyte-conditioned medium obtained both from non-irradiated (KCM-) and irradiated keratinocytes (KCM+). Again, differentially expressed genes were obtained with SAM and we performed a pathway enrichment analysis ([288]Table 6). We did not observe any significantly downregulated pathways in melanocytes growing with KCM+ vs KCM-. As regards the upregulated pathways, various pathways were affected in LM, most of them related to signaling pathways. We did not detect any upregulated pathway in DM, which suggests that DMs could have lower requirements for keratinocyte-derived factors to start the response mechanisms against UV irradiation. On the contrary, LM show a significant upregulation of several pathways when cultivated with KCM+ compared to KCM-, which could suggest that for these cells, the type or the concentration of factors present in KCM- is not enough and they require more factors to engage in certain metabolic activities. Table 6. KEGG pathway enrichment analysis for genes upregulated in the culture with KCM+ vs KCM- (significant pathways for the downregulated ones were not observed). LM DM Pathway Adj-p Pathway Adj-p 6h Ubiquitin mediated proteolysis 1.00E-03 - mTOR signaling pathway 4.34E-02 12h Neurotrophin signaling pathway 8.03E-05 - Phosphatidylinositol signaling system 5.50E-03 Endocytosis 2.28E-02 24h RNA transport 1.50E-03 Focal adhesion 4.85E-05 Insulin signaling pathway 4.70E-03 Ubiquitin mediated proteolysis 4.90E-03 Homologous recombination 7.60E-03 mRNA surveillance pathway 8.30E-03 SNARE interactions in vesicular transport 9.50E-03 Spliceosome 1.24E-02 Cell cycle 1.66E-02 Pyrimidine metabolism 1.89E-02 Ribosome biogenesis in eukaryotes 2.23E-02 Wnt signaling pathway 4.46E-02 [289]Open in a new tab Among the results obtained, of particular interest was the mTOR signaling pathway, which was upregulated in LM at 6h after UVB irradiation growing in KCM+. mTOR can be activated by UVR through the triggering of growth factor receptors bearing receptor tyrosine kinase (RTK) activity [[290]56–[291]58] like keratinocyte-derived EGF, FGF or HGF. mTOR signaling reciprocally interacts with p53 as a life/death regulator of irradiated skin cells. It has been shown that upon activation by UVR, mTOR can inhibit apoptosis and force cell cycle transition, or drive cells into senescence. This work reveals that the keratinocyte-derived factors activate the mTOR signaling pathway in LM to induce cell proliferation, consistent with the upregulation of cell cycle observed later at 24 hours. We propose that in this case mTOR forces cell cycle transition. This, however, could increase the susceptibility to develop melanoma, especially if DNA damage caused by UVB has not been repaired yet. In fact, mTOR pathway has been shown to be activated in the majority of malignant melanomas [[292]59]. The fact that this pathway was activated in LM in culture with KCM+ suggests that some keratinocyte-derived factors, secreted after the irradiation of keratinocytes with UVB, could also be at the base of melanocytes’ malignancy. Other signaling pathways that were upregulated in the presence of KCM+ are also activated by keratinocyte derived factors, such as the neurotrophin signaling pathway, which is activated by NGF and promotes the survival of melanocytes. Identification of differentially expressed genes in LM and DM under basal conditions In order to identify putative candidate genes involved in normal pigmentation variability, we compared the transcriptional profiles of DM and LM under basal conditions (i.e. at time 0, without irradiation). No significantly overrepresented pathways were observed here. Therefore, we focused on the 50 most significant genes in each category (Tables [293]7 and [294]8). The most significant genes upregulated in LM ([295]Table 7) were ATP6V0B and ATP6VOD1. These encode two components of the V-ATPase, which is responsible for maintaining an adequate acidic environment within melanosomes for the synthesis of melanin [[296]60]. The most significantly upregulated gene in DM compared to LM ([297]Table 8) was MIF. MIF has been identified as a regulator of melanogenesis, as it shows D-dopachrome tautomerase activity, which transforms D-dopachrome, dopaminechrome or its derivatives into precursors of melanin or neuromelanin [[298]61]. It has also been suggested that MIF mediates melanogenesis in the skin through the activation of protease-activated receptor (PAR-2) and stem cell factor (SCF) expression in keratinocytes after exposure to UVB [[299]62]. Interestingly, Polimanty et al [[300]63] reported a correlation between the CNV 22q11.23 containing the gene MIF with environmental variables. In particular, they suggested that MIF-related gene dosage could be associated with the adaptation to UVR, and that darker skins were correlated with haplotypes carrying no deletions. Copy number variability, and the higher frequency of deletions at this locus in light skinned individuals could be leading to a decreased MIF gene dosage, as observed in this work. Table 7. Top 50 upregulated genes in LM vs DM under basal conditions (non- coding RNAs are indicated with *) (bonferroni-adjusted p-value <0.0001). Locus name Accession number Description ATP6V0B [301]NM_004047 ATPase, H+ transporting, lysosomal 21kDa, V0 subunit b ATP6V0D1 [302]NM_004691 ATPase, H+ transporting, lysosomal 38kDa, V0 subunit d1 FUT6 [303]NM_000150 fucosyltransferase 6 (alpha (1,3) fucosyltransferase) SLC16A12 [304]NM_213606 solute carrier family 16, member 12 HSCB [305]NM_172002 HscB iron-sulfur cluster co-chaperone homolog (E. coli) ZNF865 [306]NM_001195605 zinc finger protein 865 EFS [307]NM_005864 embryonal Fyn-associated substrate KRT31 [308]NM_002277 keratin 31 RPL36A-HNRNPH2 [309]NM_001199973 RPL36A-HNRNPH2 readthrough JMJD5 [310]NM_001145348 jumonji domain containing 5 HIST1H4C [311]NM_003542 histone cluster 1, H4c GFRA3 [312]NM_001496 GDNF family receptor alpha 3 KIAA1826 [313]NM_032424 KIAA1826 DNAJC19 [314]NM_145261 DnaJ (Hsp40) homolog, subfamily C, member 19 HAP1 [315]NM_177977 huntingtin-associated protein 1 UXS1 [316]NM_025076 UDP-glucuronate decarboxylase 1 SCAMP1 [317]NM_004866 secretory carrier membrane protein 1 LTA4H [318]NM_000895 leukotriene A4 hydrolase DYNC1LI1 [319]NM_016141 dynein, cytoplasmic 1, light intermediate chain 1 LRRFIP2 [320]NM_006309 leucine rich repeat (in FLII) interacting protein 2 C10orf88 [321]NM_024942 chromosome 10 open reading frame 88 PDCD1 [322]NM_005018 programmed cell death 1 MZT2A ENST00000491265 mitotic spindle organizing protein 2A MCM3 [323]NM_002388 minichromosome maintenance complex component 3 C6orf163 [324]NM_001010868 chromosome 6 open reading frame 163 DTX4 [325]NM_015177 deltex homolog 4 (Drosophila) CLCN6 [326]NM_021735 chloride channel 6 (CLCN6) RFK [327]NM_018339 riboflavin kinase WHSC2 [328]NM_005663 Wolf-Hirschhorn syndrome candidate 2 FGFRL1 [329]NM_001004356 fibroblast growth factor receptor-like 1 BTBD6 [330]NM_033271 BTB (POZ) domain containing 6 N4BP1 [331]NM_153029 NEDD4 binding protein 1 MAP2K6 [332]NM_002758 mitogen-activated protein kinase kinase 6 POMP [333]NM_015932 proteasome maturation protein GABPA [334]NM_002040 GA binding protein transcription factor, alpha subunit 60kDa UPF3A [335]NM_023011 UPF3 regulator of nonsense transcripts homolog A (yeast) PLEKHA3 [336]NM_019091 pleckstrin homology domain containing, family A member 3 CD276 [337]NM_001024736 CD276 molecule (CD276) ENTPD2 [338]NM_203468 ectonucleoside triphosphate diphosphohydrolase 2 DEDD [339]NM_032998 death effector domain containing FAM70B ENST00000375348 family with sequence similarity 70, member B MCM5 [340]NM_006739 minichromosome maintenance complex component 5 LOC100131257* [341]NR_034022 zinc finger protein 655 pseudogene SCARNA13 [342]NR_003002 small Cajal body-specific RNA 13 SMNDC1 [343]NM_005871 survival motor neuron domain containing 1 CALML4 [344]NM_033429 calmodulin-like 4 C1orf131 [345]NM_152379 chromosome 1 open reading frame 131 RNGTT [346]NM_003800 RNA guanylyltransferase and 5'-phosphatase KCNQ3 [347]NM_004519 potassium voltage-gated channel, KQT-like subfamily, member 3 WASF3 [348]NM_006646 WAS protein family, member 3 [349]Open in a new tab Table 8. Top 50 upregulated genes in DM vs LM under basal conditions (non- coding RNAs are indicated with *) (bonferroni-adjusted p-value <0.0001). Locus name Accession number Description MIF [350]NM_002415 macrophage migration inhibitory factor TTC19* ENST00000395886 tetratricopeptide repeat domain 19 CYTB ENST00000361789 mitochondrially encoded cytochrome b NBEA [351]NM_015678 neurobeachin CTSO [352]NM_001334 cathepsin O SNORA23 [353]NR_002962 small nucleolar RNA, H/ACA box 23 PMP22 [354]NM_000304 peripheral myelin protein 22 CXCL1 [355]NM_001511 chemokine ligand (melanoma growth stimulating activity, alpha) CDKN2A [356]NM_058197 cyclin-dependent kinase inhibitor 2A (melanoma, p16) LDB3 [357]NM_001171610 LIM domain binding 3 MIPEP [358]NM_005932 mitochondrial intermediate peptidase MAPK8 [359]NM_139047 mitogen-activated protein kinase 8 SNORD15A [360]NR_000005 small nucleolar RNA, C/D box 15A SNAR-A3* [361]NR_024214 small ILF3/NF90-associated RNA A3 ASCC1 [362]NM_015947 activating signal cointegrator 1 complex subunit 1 ZNF235 [363]NM_004234 zinc finger protein 235 MBIP [364]NM_001144891 MAP3K12 binding inhibitory protein 1 C13orf38 [365]NM_001198908 chromosome 13 open reading frame 38 LOC100132707* [366]NR_024477 hypothetical LOC100132707 UTRN [367]NM_007124 utrophin CALM2 [368]NM_001743 calmodulin 2 (phosphorylase kinase, delta) MOCS1* [369]NM_005943 molybdenum cofactor synthesis 1 ZNF212 [370]NM_012256 zinc finger protein 212 KIAA0090 [371]NM_015047 KIAA0090 (KIAA0090) SNORD3B-1 [372]NR_003271 small nucleolar RNA, C/D box 3B-1 HLX [373]NM_021958 H2.0-like homeobox C9orf72 [374]NM_145005 chromosome 9 open reading frame 72 SEC23B [375]NM_032985 Sec23 homolog B (S. cerevisiae) WRAP73 [376]NM_017818 WD repeat containing, antisense to TP73 MAN1A1 [377]NM_005907 mannosidase, alpha, class 1A, member 1 S100B [378]NM_006272 S100 calcium binding protein B CCDC93 [379]NM_019044 coiled-coil domain containing 93 ZNF3 [380]NM_032924 zinc finger protein 3 FTH1 [381]NM_002032 ferritin, heavy polypeptide 1 RAB30 [382]NM_014488 RAB30, member RAS oncogene family RDM1 [383]NM_001034836 RAD52 motif 1 BTAF1 [384]NM_003972 BTAF1 RNA polymerase II, HLA-F [385]NM_018950 major histocompatibility complex, class I, F CABYR [386]NM_012189 calcium binding tyrosine-(Y)-phosphorylation regulated MAP4K2 [387]NM_004579 mitogen-activated protein kinase 2 PRPF18 ENST00000320054 PRP18 pre-mRNA processing factor 18 homolog (S. cerevisiae) CALM3 [388]NM_005184 calmodulin 3 (phosphorylase kinase, delta) ALKBH4 [389]NM_017621 alkB, alkylation repair homolog 4 (E. coli) LOC399744* [390]NR_024497 hypothetical LOC399744 SETD6 [391]NM_024860 SET domain containing 6 SH3TC2 [392]NM_024577 SH3 domain and tetratricopeptide repeats 2 ARHGAP35 [393]NM_004491 Rho GTPase activating protein 35 CHCHD6 [394]NM_032343 coiled-coil-helix-coiled-coil-helix domain containing 6 RC3H2 [395]NM_018835 ring finger and CCCH-type domains 2 WDR46 [396]NM_005452 WD repeat domain 46 [397]Open in a new tab For the other genes in Tables [398]7 and [399]8 we did not find any evident direct correlation with pigmentary phenotype. Validation by RT-qPCR Six genes were selected for validation of the microarrays' results, which showed either a change of expression after UV treatment or a differential expression between LM and DM (ATP6VOB, TP53I3, MDM2, MIF, RPL6 and FDXR), and measured their expression levels by quantitative real-time PCR (RT-qPCR). We assessed the expression of 4 melanocytic cell lines (2 DM and 2 LM) at basal conditions and at 6 and 12 hours after UVB irradiation. The expression patterns and direction of changes of all of the genes were consistent with the microarray data ([400]Fig 4), observing a significant increase in the expression of TP53I3, MDM2, RPL6 and FDXR in both LM and DM after UVB. The expression analysis of ATP6VOB and MIF also supported the differential expression of these genes by LM and DM, being ATP6VOB more expressed by LM, while MIF was more significantly expressed by DM (both at basal conditions and after UVB irradiation). Fig 4. Gene expression of genes FDX6, RPL6, MDM2, TP53I3, ATP6VOB and MIF assessed by RT-qPCR. [401]Fig 4 [402]Open in a new tab Unpaired t-test; *** p<0.0001; ** p<0.001; * p<0.01. Natural selection tests In order to assess the biological relevance of the genes that were differentially expressed between DM and LM under basal conditions, we performed evolutionary neutrality tests on these genes (Tables [403]7 and [404]8) using the populations from the 1000 Genomes Project (1KGP). For this, we performed a first screening of different neutrality tests using the 1000 Genomes Selection Browser to identify putative signatures of selection. After multiple test correction, the gene ATP6V0D1 (ATPase, H+ Transporting, Lysosomal 38kDa, V0 Subunit D1) seemed to deviate from neutrality in the European populations. Further neutrality tests using DnaSP [[405]28] supported significant signatures of selection acting on ATP6V0D1 in Europeans (Tajima's D: -2.31, p-value = 0; Fay & Wu's H: -10.66, p-value = 0.001), thus suggesting that this gene might be involved in human pigmentary phenotype. This reinforces the notion that selective pressures can shape pigmentation variability by driving the evolution of melanosomal genes. So, besides the well-known OCA2, SLC45A2 and SLC24A5, we support ATP6V0D1 as an additional melanosomal-membrane gene that has been subjected to selective pressures and might be involved in pigmentation variability in Europeans. No deviations from neutrality were detected in any population for the MIF gene (data not shown). However, we should take into account that MIF is embedded in a CNV [[406]63] and in a previous work we observed how a variation in copy number can interfere with neutrality tests by altering the frequencies of polymorphisms leading to an excess of detected homozygosity [[407]64]. A loss of copies would result in apparent homozygosity, and duplications of one allele would mask possible variant alleles in sequencing or genotyping experiments. Therefore, although with the available tools and our knowledge we cannot detect deviations from neutrality, we cannot still exclude the possibility that this gene is under selection. Conclusions We have provided an overview of the most significant genes that are up and downregulated in response to UVB irradiation and revealed the interaction of ribosomal proteins and P53 signaling pathway in the response to UVB in both DM and LM. We have also observed that DM and LM show differentially expressed genes after irradiation and in particular in the first 6 hours. These are mainly associated with inflammatory skin reactions, cell survival or melanoma. Furthermore, the culture with KCM+ compared with KCM- had a noticeable effect on LM, but not in DM, triggering various signaling pathways in LM such as the mTOR signaling pathway. And importantly, the comparison of the transcriptional profile of LM and DM under basal conditions allowed us to highlight the significant involvement of MIF and ATP6V0B in the normal variability of human skin pigmentation. Supporting Information S1 Table. Top 50 upregulated genes in LM vs DM 12 hours after UVB. (PDF) [408]Click here for additional data file.^ (124.8KB, pdf) S2 Table. Top 50 upregulated genes in DM vs LM 12 hours after UVB. (PDF) [409]Click here for additional data file.^ (123.7KB, pdf) S3 Table. Top 50 upregulated genes in LM vs DM 24 hours after UVB. (PDF) [410]Click here for additional data file.^ (124.6KB, pdf) S4 Table. Top 50 upregulated genes in DM vs LM 24 hours after UVB. (PDF) [411]Click here for additional data file.^ (48.9KB, pdf) Acknowledgments