Figures
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.
Citation: López S, Smith-Zubiaga I, García de Galdeano A, Boyano MD, García O, Gardeazábal J, et al. (2015) Comparison of the Transcriptional Profiles of Melanocytes from Dark and Light Skinned Individuals under Basal Conditions and Following Ultraviolet-B Irradiation. PLoS ONE 10(8): e0134911. https://doi.org/10.1371/journal.pone.0134911
Editor: Thomas G. Hofmann, German Cancer Research Center, GERMANY
Received: June 2, 2015; Accepted: July 16, 2015; Published: August 5, 2015
Copyright: © 2015 López et al. This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited
Data Availability: All microarray data are available from the Gene Expression Omnibus database (Accession Number GSE70280).
Funding: This work was supported by the former Spanish Ministerio de Ciencia e Innovación, project CGL2008-04066/BOS to S.A.; by the Dept. Educación, Universidades e Investigación of the Basque Government, project IT542-10 to C.R.; the University of the Basque Country program UFI11/09; a predoctoral fellowship from the Dept. Educación, Universidades e Investigación of the Basque Government to S.L. (BFI09.248). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.
Competing interests: The authors have declared that no competing interests exist.
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 [1–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 [5–7]. One of the first reports using cDNA microarrays in various cell lines of human melanocytes for around 9,000 human genes [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/cm2). Later, Yang et al. [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. [5], Yang et al. [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 [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 [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 [10–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+) [19–20], in other studies the medium is collected from keratinocytes that have not been previously irradiated (hereinafter KCM-) [10, 21–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% CO2 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/cm2 UVB, based on our previous work [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/cm2 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 (Fig 1).
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 [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 [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, [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 [26], adjusted for multiple tests.
Pathway enrichment analysis
Enrichment analysis was performed using Web-based Gene Set Analysis Toolkit (WebGestalt) (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 r2>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 (http://hsb.upf.edu/) [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 [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) [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 [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 [30], that has been used for the analysis of microarray data in search of outlier genes [31] or to identify temporal patterns in time-series analyses [32]. The PCA (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.
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 [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 1 and 2). The adjusted p-value for all these genes was <0.0001.
Common upregulated genes after UVB irradiation.
Some of the genes included in this category (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 [33] and melanocytes [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 [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. [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 (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 [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 [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 [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 (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 S1–S4 Tables).
Upregulated genes in LM vs DM 6 hours after UVB irradiation.
The significant upregulation of EDA2R in LM (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 [38–39], and its paralog, EDAR, has been strongly associated with skin pigmentation variability in humans [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 [41].
We also identified as upregulated some genes related to melanoma, like CDKN2A, whose expression has already been reported to be induced by UV radiation [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 [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 (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 [44], PTGS2 [45] or CCL2 [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.
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 [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 (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).
The role of P53, a tumour suppressor that promotes either cell cycle arrest and DNA repair, apoptosis or senescence [48] in the response to UVB has already been reported [6]. Our results are consistent with the proposed mechanism of P53 pathway regulation by ribosomal proteins [49–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 (Fig 3). Among the ribosomal proteins that can bind to MDM2 are RPL5 [52] and RPL11 [53], both of which were among the upregulated ribosomal genes in this work.
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 (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 [54], thus preventing the development and progression of melanoma [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 (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.
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 [56–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 [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 7 and 8). The most significant genes upregulated in LM (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 [60]. The most significantly upregulated gene in DM compared to LM (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 [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 [62]. Interestingly, Polimanty et al [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.
For the other genes in Tables 7 and 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 (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).
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 7 and 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 [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 [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 [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.
https://doi.org/10.1371/journal.pone.0134911.s001
(PDF)
S2 Table. Top 50 upregulated genes in DM vs LM 12 hours after UVB.
https://doi.org/10.1371/journal.pone.0134911.s002
(PDF)
S3 Table. Top 50 upregulated genes in LM vs DM 24 hours after UVB.
https://doi.org/10.1371/journal.pone.0134911.s003
(PDF)
S4 Table. Top 50 upregulated genes in DM vs LM 24 hours after UVB.
https://doi.org/10.1371/journal.pone.0134911.s004
(PDF)
Acknowledgments
The authors thank for technical and human support provided by the Genomics and Proteomics Service of the UPV/EHU (SGIker).
Author Contributions
Conceived and designed the experiments: SL ISZ SA. Performed the experiments: SL ISZ AGG SA. Analyzed the data: SL ISZ SA. Contributed reagents/materials/analysis tools: MDB OG JG CMC NI CR AGG. Wrote the paper: SL ISZ SA.
References
- 1. Young AR. The sunburn cell. Photodermatol. 1987;4: 127–134. pmid:3317295
- 2. Kricker A, Armstrong BK, English DR. Sun exposure and non-melanocytic skin cancer. Cancer Cause Control. 1994;5: 367–392.
- 3. Elwood JM. Melanoma and sun exposure. Semin Oncol. 1996;23: 650–666. pmid:8970584
- 4. Povey JE, Darakhshan F, Robertson K, Bisset Y, Mekky M, Rees J, et al. DNA repair gene polymorphisms and genetic predisposition to cutaneous melanoma. Carcinogenesis. 2007;28: 1087–1093. pmid:17210993
- 5. Valéry C, Grob JJ, Verrando P. Identification by cDNA microarray technology of genes modulated by artificial ultraviolet radiation in normal human melanocytes: relation to melanocarcinogenesis. J Invest Dermatol. 2001;117: 1471–1482. pmid:11886511
- 6. Yang G, Zhang G, Pittelkow MR, Ramoni M, Tsao H. Expression profiling of UVB response in melanocytes identifies a set of p53-target genes. J Invest Dermatol. 2006;126: 2490–2506. pmid:16888633
- 7. Alonso S, Izagirre N, Smith-Zubiaga I, Gardeazabal J, Díaz-Ramón JL, Díaz-Pérez JL, et al. Complex signatures of selection for the melanogenic loci TYR, TYRP1 and DCT in humans. BMC Evol Biol. 2008;8: 74. pmid:18312627
- 8. Tadokoro T, Yamaguchi Y, Batzer J, Coelho SG, Zmudzka BZ, Miller SA, et al. Mechanisms of skin tanning in different racial/ethnic groups in response to ultraviolet radiation. J Invest Dermatol. 2005;124: 1326–1332. pmid:15955111
- 9. Haltaufderhyde KD, Oancea E. Genome-wide transcriptome analysis of human epidermal melanocytes. Genomics. 2014;104: 482–489. pmid:25451175
- 10. Gordon PR, Mansur CP, Gilchrest BA. Regulation of human melanocyte growth, dendricity, and melanization by keratinocyte derived factors. J Invest Dermatol. 1989;92: 565–572. pmid:2467948
- 11. Yaar M, Grossman K, Eller M, Gilchrest BA. Evidence for nerve growth factor-mediated paracrine effects in human epidermis. J Cell Biol. 1991;115: 821–828. pmid:1655813
- 12. Chakraborty AK, Funasaka Y, Slominski A, Ermak G, Hwang J, Pawelek JM, et al. Production and release of proopiomelanocortin (POMC) derived peptides by human melanocytes and keratinocytes in culture: regulation by ultraviolet B. Biochim Biophys Acta. 1996;1313: 130–138. pmid:8781560
- 13. Hirobe T. Endothelins are involved in regulating the proliferation and differentiation of mouse epidermal melanocytes in serum-free primary culture. J Invest Dermatol Symp Proc. 2001;6: 25–31.
- 14. Hirobe T, Osawa M, Nishikawa S. Steel factor controls the proliferation and differentiation of neonatal mouse epidermal melanocytes in culture. Pigm Cell Res. 2003;16: 644–655.
- 15. Hirobe T, Osawa M, Nishikawa S. Hepatocyte growth factor controls the proliferation of cultured epidermal melanoblasts and melanocytes from newborn mice. Pigm Cell Res. 2004;17: 51–61.
- 16. Scott G, Leopardi S, Printup S, Malhi N, Seiberg M, Lapoint R. Proteinase-activated receptor-2 stimulates prostaglandin production in keratinocytes: analysis of prostaglan- din receptors on human melanocytes and effects of PGE2 and PGF2a on melanocyte dendricity. J Invest Dermatol. 2004;122: 1214–1224. pmid:15140225
- 17. Yamaguchi Y, Brenner M, Hearing VJ. The regulation of skin pigmentation. J Biol Chem. 2007;282: 27557–27561. pmid:17635904
- 18. Abdel-Malek ZA, Kadekaro AL, Swope VB. Stepping up melanocytes to the challenge of UV exposure. Pigm Cell Melanoma R. 2010;23: 171–186.
- 19. Jamal S, Schneider RJ. UV-induction of keratinocyte endothelin-1 downregulates E-cadherin in melanocytes and melanoma cells. J Clin Invest. 2002;110: 443–452. pmid:12189238
- 20. López S, Alonso S, García de Galdeano A, Smith-Zubiaga I. Melanocytes from dark and light skin respond differently after ultraviolet B irradiation: effect of keratinocyte conditioned medium. Photodermatol Photoimmunol Photomed. 2015;31: 149–158. pmid:25740555
- 21. Abdel-Naser MB. Differential effects on melanocyte growth and melanization of low vs. high calcium keratinocyte-conditioned medium. Br J Dermatol. 1999;140: 50–55. pmid:10215767
- 22. Maeda K, Tomita Y. Mechanism of the Inhibitory Effect of Tranexamic Acid on Melanogenesis in Cultured Human Melanocytes in the Presence of Keratinocyte-conditioned Medium. J Health Sci. 2007;53: 389–396.
- 23. Motakis ES, Nason GP, Fryzlewicz P, Rutter GA. Variance stabilization and normalization for one-color microarray data using a data-driven multiscale approach. Bioinformatics. 2006;22: 2547–2453. pmid:16877753
- 24. Bolstad BM, Irizarry RA, Astrand M, Speed TP. A Comparison of Normalization Methods for High Density Oligonucleotide Array Data Based on Bias and Variance. Bioinformatics. 2003;19: 185–193. pmid:12538238
- 25. Tusher VG, Tibshirani R, Chu G. Significance analysis of microarrays applied to the ionizing radiation response. Proc Natl Acad Sci U S A. 2001;98: 5116–5121. pmid:11309499
- 26. Storey JD. A direct approach to false discovery rates. J Roy Stat Soc B. 2002;64: 479–498.
- 27. Pybus M, Dall'Olio GM, Luisi P, Uzkudun M, Carreño-Torres A, Pavlidis P, et al. 1000 Genomes Selection Browser 1.0: a genome browser dedicated to signatures of natural selection in modern humans. Nucleic Acids Res. 2014;42: D903–909. pmid:24275494
- 28. Librado P, Rozas J. DnaSP v5: a software for comprehensive analysis of DNA polymorphism data. Bioinformatics. 2009;25: 1451–1452. pmid:19346325
- 29. Amigo J, Salas A, Phillips C, Carracedo A. SPSmart: adapting population based SNP genotype databases for fast and comprehensive web access. BMC Bioinformatics. 2008;9: 428. pmid:18847484
- 30.
Basilevsky A. Statistical Factor Analysis and Related Methods, Theory and Applications. John Wiley & Sons. New York; 1994.
- 31. Hilsenbeck SG, Friedrichs WE, Schiff R, O’Connell P, Hansen RK, Osborne CK, et al. Statistical Analysis of Array Expression Data as Applied to the Problem of Tamoxifen Resistance. J Natl Cancer Institute. 1999;91: 453–459.
- 32. Chu S, DeRisi J, Eisen M, Mulholland J, Botstein D, Brown PO, et al. The transcriptional program of sporulation in budding yeast. Science. 1998;282: 699–705. pmid:9784122
- 33. de la Fuente H, Lamana A, Mittelbrunn M, Perez-Gala S, Gonzalez S, García-Diez A, et al. Identification of genes responsive to solar simulated UV radiation in human monocyte-derived dendritic cells. PLoS ONE. 2009;4: e6735. pmid:19707549
- 34. Fukunaga-Kalabis M, Martinez G, Liu ZJ, Kalabis J, Mrass P, Weninger W, et al. CCN3 controls 3D spatial localization of melanocytes in the human skin through DDR1. J Cell Biol. 2006;175: 563–569. pmid:17101694
- 35. Braeuer RR, Zigler M, Kamiya T, Dobroff AS, Huang L, Choi W, et al. Galectin-3 Contributes to Melanoma Growth and Metastasis via Regulation of NFAT1 and Autotaxin. Cancer Res. 2012;15: 5757–5766.
- 36. Baumann L. How to Prevent Photoaging? J Invest Dermatol. 2005;125: 7–13.
- 37. Yaar M, Park HY. Melanocytes: a window into the nervous system. J Invest Dermatol. 2012;132: 835–845. pmid:22158549
- 38. Sabeti PC, Varilly P, Fry B, Lohmueller J, Hostetter E, Cotsapas C, et al. Genome-wide detection and characterization of positive selection in human populations. Nature. 2007;18: 913–918.
- 39. Hilmer AM, Freudenberg J, Myles S, Herms S, Tang K, Hughes DA, et al. Recent positive selection of a human androgen receptor/ectodysplasin A2 receptor haplotype and its relationship to male pattern baldness. Hum Genet. 2009;126: 255–264. pmid:19373488
- 40. Hider JL, Gittelman RM, Shah T, Edwards M, Rosenbloom A, Akey JM, et al. Exploring signatures of positive selection in pigmentation candidate genes in populations of East Asian ancestry. BMC Evol Biol. 2013;13: 150. pmid:23848512
- 41. Casto AM, Henn BM, Kidd JM, Bustamante CD, Feldman MW. A tale of two haplotypes: the EDA2R/AR Intergenic region is the most divergent genomic segment between Africans and East Asians in the human genome. Hum Biol. 2012;84: 641–694. pmid:23959643
- 42. Parvey S, Conroy S, Russell T, Gabrielly B. Ultraviolet radiation induces p16 expression in human skin. Cancer Res. 1999;59: 4185–4189. pmid:10485451
- 43. Rogers SL, Karcher RL, Roland JT, Minin AA, Steffen W, Gelfand VI. Regulation of melanosome movement in the cell cycle by reversible association with myosin V. J Cell Biol. 1999;146: 1265–1276. pmid:10491390
- 44. Urbanski A, Schwarz T, Neuner P, Krutmann J, Kimbauer R, Köck A, et al. Ultraviolet light induces increased circulating interleukin-6 in humans. J Invest Dermatol. 1990;94: 808–811. pmid:2355183
- 45. Buckman SY, Gresham A, Hale P, Hruza G, Anast J, Masferrer J, et al. COX-2 expression is induced by UVB exposure in human skin: implications for the development of skin cancer. Carcinogenesis. 1998;19: 723–729. pmid:9635856
- 46. Purwar R, Wittmann M, Zwirner J, Oppermann M, Kracht M, Dittrich-Breiholz O, et al. Induction of C3 and CCL2 by C3a in keratinocytes: a novel autocrine amplification loop of inflammatory skin reactions. J Immunol. 2006;177: 4444–4450. pmid:16982879
- 47. Zhou Y, Zhong Y, Wang Y, Zhang X, Batista DL, Gejman R, et al. Activation of p53 by MEG3 non-coding RNA. J Biol Chem. 2007;282: 24731–24742. pmid:17569660
- 48. Rufini A, Tucci P, Celardo I, Melino G. Senescence and aging: the critical roles of p53. Oncogene. 2013;32: 5129–5143. pmid:23416979
- 49. Zhang Y, Lu H. Signaling to p53: ribosomal proteins find their way. Cancer Cell. 2009;16: 369–377. pmid:19878869
- 50. Deisenroth C, Zhang Y. Ribosome biogenesis surveillance: probing the ribosomal protein-Mdm2-p53 pathway. Oncogene. 2010;29: 4253–4260. pmid:20498634
- 51. Donati G, Montanaro L, Derenzini M. Ribosome biogenesis and control of cell proliferation: p53 is not alone. Cancer Res. 2012;72: 1602–1607. pmid:22282659
- 52. Dai MS, Lu H. Inhibition of MDM2-mediated p53 ubiquitination and degradation by ribosomal protein L5. J Biol Chem. 2004;279: 44475–44482. pmid:15308643
- 53. Lohrum MA, Ludwig RL, Kubbutat MH, Hanlon M, Vousden KH. Regulation of HDM2 activity by the ribosomal protein L11. Cancer Cell. 2003;3: 577–587. pmid:12842086
- 54. Garrod DR. Desmosomes and hemidesmosomes. Curr Opin Cell Biol. 1993;5: 30–40. pmid:8448028
- 55. Haass NK, Herlyn M. Normal human melanocyte homeostasis as a paradigm for understanding melanoma. J Investig Dermatol Symp Pro. 2005;10: 153–163.
- 56. Jean C, Hernandez-Pigeon H, Blanc A, Charveron M, Laurent G. Epidermal growth factor receptor pathway mitigates UVA-induced G2/M arrest in keratinocyte cells. J Invest Dermatol. 2007;127: 2418–2424. pmid:17495959
- 57. Cao C, Lu S, Kivlin R, Wallin B, Card E, Bagdasarian A, et al. AMP-activated protein kinase contributes to UV- and H2O2-induced apoptosis in human skin keratinocytes. J Biol Chem. 2008;283: 28897–28908. pmid:18715874
- 58. He YY, Council SE, Feng L, Chignell CF. UVA-induced cell cycle progression is mediated by a disintegrin and metalloprotease/epidermal growth factor receptor/AKT/Cyclin D1 pathways in keratinocytes. Cancer Res. 2008;68: 3752–3758. pmid:18483258
- 59. Karbowniczek M, Spittle CS, Morrison T, Wu H, Henske EP. mTOR is activated in the majority of malignant melanomas. J Invest Dermatol. 2008;128: 980–987. pmid:17914450
- 60. Dooley CM, Schwarz H, Mueller K, Mongera A, Konantz M, Neuhauss SC, et al. Slc45a2 and V-ATPase are regulators of melanosomal pH homeostasis in zebrafish, providing a mechanism for human pigment evolution and disease. Pigment Cell Melanoma Res. 2013;26: 205–217. pmid:23205854
- 61. Matsunaga J, Sinha D, Solano F, Santis C, Wistow G, Hearing V. Macrophage migration inhibitory factor (MIF)—its role in catecholamine metabolism. Cell Mol Biol. 1999;45: 1035–1040. pmid:10644007
- 62. Enomoto A, Yoshihisa Y, Yamakoshi T, Ur Rehman M, Norisugi O, Hara H, et al. UV-B radiation induces macrophage micgration inhibitory factor-mediated melanogenesis through activation of protease-activated receptor-2 and stem cell factor in keratinocytes. Am J Pathol. 2011;178: 679–687. pmid:21281800
- 63. Polimanti R, Piacentini S, Iorio A, De Angelis F, Kozlov A, Novelletto A, et al. Haplotype differences for copy number variants in the 22q11.23 region among human populations: a pigmentation-based model for selective pressure. Eur J Hum Genet. 2015;23: 116–123. pmid:24667780
- 64. López S, García I, Smith I, Sevilla A, Izagirre N, de la Rúa C, et al. Discovery of Copy Number Variants by Multiplex Amplifiable Probe Hybridization (MAPH) in candidate pigmentation genes. Ann Hum Biol. 2014;24: 1–9.