Potential molecular characteristics in situ in response to repetitive UVB irradiation.

BACKGROUND: To identify molecular characteristics in situ in response to repetitive UVB (ultraviolet-B) irradiation.
METHODS: Microarray data from the Gene Expression Omnibus were re-analyzed to identify DEGs (differentially expressed genes) between UVB-irradiated and non-irradiated skin biopsies. Enrichment and annotation analyses were performed respectively using DAVID, and TSGene and TAG databases. PPIs (protein-protein interactions) were analyzed using STRING, and miRNAs (microRNAs) and TFs (transcription factors) were predicted separately by miRNA-related databases and ENCODE. Accordingly, the PPI network and regulatory networks were visualized using Cytoscape, and they were merged together to obtain an integrated network for mining densely connected modules.
RESULTS: Altogether, 151 up- and 64 down-regulated genes were identified between UVB-irradiated and non-irradiated skin biopsies, among which down-regulated DNAJB4 and SLIT2 were annotated as tumor-suppressors and up-regulated KIT was annotated as an oncogene. The up-regulated DEGs were significantly enriched in biological processes related to pigmentation (DCT, SOX10, TYRP1, TYR, MLPH, KIT and GPR143), while the down-regulated DEGs were dramatically related to haemopoiesis and the immune system (GPR183, INHBA, PTPRC, PLEK, CD8A and IKZF1). Furthermore, many miRNAs were screened for the DEGs, including miR-206 and miR-496 targeting KIT, miR-184 targeting DCT, and highly significant miR-337-5p, miR-21 and miR-16. Additionally, TFs were identified for the DEGs, among which PAX5 and HNF4A targeted MLPH and GPR143, respectively, while BATF, SPI1 and EP300 jointly target GPR183, PTPRC and PLEK.
CONCLUSIONS: The pigmentation and immune system implicated by DEGs, miRNAs and TFs might be important molecular mechanisms in response to UVB irradiation.

Background

Sunlight is essential for human life, but a previous study has revealed that it may induce 2 million new cases of non-melanoma skin cancer every year in the United States [1]. UV (ultraviolet) may impair the ability of basal keratinocytes to remove melanin and thus lead to hyperpigmentation [2]. Furthermore, UV radiation might be the major contributor to the risk of certain types of skin cancers in UV-radiated outdoor workers, such as basal cell carcinoma and squamous cell carcinoma [3]. UVB (ultraviolet-B) is a type of UV with wavelength ranging from 280 to 320 nm, and excessive UVB irradiation from sunlight may induce oxidative damage and inflammation in skin that further causes sunburn, photo-aging and various skin cancers [4].

Numerous researchers have been exploring the molecular characteristics in response to UVB exposure using animal models or cell culture systems. The activity of tyrosinase essential for melanogenic cascade can be stimulated by UVB exposure that further increases melanin synthesis required for UVB-induced delayed tanning [5]. The membrane-bound stem cell factor is up-regulated and then activates neighboring melanocytes through kit receptors in the course of UVB-caused pigmentation [6]. Nrf2 (nuclear factor, erythroid 2-like 2) may exert protective effects on an inflammation and a sunburn reaction that were induced by UVB exposure via inhibiting inflammation and extracellular matrix degradation [7]. More dangerously, UVB exposure can down-regulate the tumor suppressor TGFβ (Transforming growth factor β) signaling, suggesting its oncogenic effects [8]. Additionally, miRNAs (microRNAs) and TFs (transcription factors) are also implicated in the UVB irradiation-induced cellular response via regulating target genes [9, 10].

Although those previous studies have facilitated our understanding of molecular mechanisms in response to UVB exposure, they might not precisely reveal the molecular characteristics in situ where human skin is repetitively irradiated by UVB. Wonseon Choi et al. have demonstrated that UVB irradiation in situ can induce a large number of significant UVB-responsive genes implicated in regulating skin pigmentation, which has a different molecular characterization in response to UVA (ultraviolet-A) exposure [11]. However valuable, the above research focuses on the different molecular characteristics of human skin pigmentation under various types of UV irradiation. It is well known that different analysis procedures can obtain different findings. From a distinct perspective, this study was designed to explore the underlying molecular mechanisms in situ in response to UVB irradiation using various bioinformatics methods.

This study re-analyzed the published microarray data [11] to screen out DEGs (differentially expressed genes) between UVB-irradiated and non-irradiated skin biopsies, which were then subjected to enrichment analyses to identify biological metabolism and annotation analysis to identify tumor-associated DEGs. Then, PPIs (protein-protein interactions) among the DEGs were screened out to construct a PPI network, whilst TFs and miRNAs were screened out to construct regulatory networks. Accordingly, an integrated network among TFs, miRNAs and DEGs was constructed and further mined to identify densely connected modules.

Methods

This article does not contain any studies with human participants or animals performed by any of the authors.

Gene expression profiles

To obtain a global view of the gene expression patterns induced by exposure to UVB, gene expression profiling data [11] from UVB-irradiated skin biopsies and matched controls were downloaded from the GEO (Gene Expression Omnibus, http://www.ncbi.nlm.nih.gov/geo/) database [12]. The skin biopsies were obtained from six volunteers with Fitzpatrick skin type II–III, whose average age was 37.3 ± 15.7 years. All donors provided written informed consent before enrollment in the irradiation experiment, which was conducted in accordance with the Helsinki guidelines and approved by the Research Involving Human Subjects Committees of Beiersdorf AG. For UVB irradiation, the backs of all participants were irradiated under 290–320 nm for two weeks with a frequency of five times per week using a custom-made filter combination (Tafelmayer, Rosenheim, Germany). Then, UVB-irradiated skin biopsies were obtained 3 days after the last irradiation, and corresponding non-irradiated skin biopsies were taken from comparable skin regions as controls. Then, total RNA was extracted from all skin biopsies to perform single-colour hybridization on an Agilent-014850 Whole Human Genome Microarray 4x44K G4112F.

Raw data preprocessing and differential expression analysis

The downloaded raw expression profiles were preprocessed through background correction, quantile normalization, probe summarization and transformation of probe IDs to gene symbols by the Limma (Linear Models for Microarray data, version 1.8.18, http://www.bioconductor.org/packages/release/bioc/html/limma.html) package in R [13]. Accordingly, a gene expression matrix was obtained and subsequently subjected to differential expression analysis between UVB-irradiated and non-irradiated skin biopsies using the t-test. The p-value of each gene was adjusted by the Benjamini & Hochberg method [14]. Genes with a P-value < 0.05 and a |log2 fold change (FC)| >0.5 were selected as DEGs. Additionally, in consideration of the complex regulatory effects of TFs and miRNAs on target genes, genes with a P-value < 0.01 were identified as expressed genes for subsequent analysis of TFs and miRNAs.

Enrichment analyses and functional annotations

Functional and pathway enrichment analyses can be implemented using the web-accessible DAVID (Database for Annotation, Visualization and Integrated Discovery, http://david.abcc.ncifcrf.gov/) v6.7 [15], which augments the biological meaning of genome-scale datasets by providing GO (Gene Ontology, version 1.2, http://www.geneontology.org/) [16] functional and KEGG (Kyoto Encyclopedia of Genes and Genomes, version 58, http://www.genome.jp/kegg/pathway.html) [17] pathway enrichment analyses. To interpret the biological meanings of DEGs, GO and KEGG enrichment analyses were respectively implemented for the significantly up-regulated and down-regulated genes, respectively, with all genes in the human genome being used as background. The thresholds were P-values < 0.05 and the numbers of enriched genes in each term >2. Additionally, tumor-associated genes were screened separately from the DEGs using TSGene (Tumor Suppressor Gene, http://bioinfo.mc.vanderbilt.edu/TSGene/) (20 March, 2015) [18] and TAG (Tumor Associated Gene, http://www.binfo.ncku.edu.tw/TAG/GeneDoc.php) (20 March, 2015) databases [19].

Protein-protein interaction network construction

To further explore interactions among the DEGs, they were inputted into the online database STRING v9.1 (Search Tool for the Retrieval of Interacting Genes, http://string-db.org/) [20]. Subsequently, interaction pairs were screened out based on a combined score threshold >0.4, and they were visualized by constructing a PPI network in Cytoscape v2.6 (http://cytoscape.org/) [21]. The hub nodes, corresponding to DEGs with predominant roles in biological processes [22], were further screened from the PPI network.

Regulatory network construction

To investigate the up-stream regulatory mechanisms, a total of 6 databases were used to predict miRNAs for the expressed genes with the thresholds of P-values < 0.05 and the number of target genes ≥2, including miRanda v1.9, MirTarget2 v2.0, PicTar v4, PITA v6, TargetScan v4.1 and miRecords v3 [23]. Moreover, TFs were screened out among the expressed genes based on the ENCODE (Encyclopedia of DNA Elements, https://www.encodeproject.org/) project consortium 2007 [24] and subsequently subjected to detection of TF activity using the partial least square algorithm [25]. Accordingly, regulatory networks of miRNAs-target genes and TFs-target genes were constructed using Cytoscape v2.6.

Integrated network construction and module screening

The PPI network and the two regulatory networks were further merged together to obtain one integrated network among TFs, miRNAs and DEGs. Subsequently, the integrated network was mined to screen out densely connected modules with default cutoffs (degree cutoff: 2, node score cutoff: 0.2, K-core: 2, and maximum depth: 100) using MCODE v1.31 (Molecular Complex Detection, http://baderlab.org/Software/MCODE) [26]. Functional enrichment analyses were performed for the densely connected modules to identify their biological relevance using DAVID v6.7.

Results

Differentially expressed genes

With cutoffs of P-value < 0.05 and |log2 FC| >0.5, altogether 215 DEGs including 151 up- and 64 down-regulated genes were identified in UVB-irradiated skin biopsies compared with non-irradiated controls. Moreover, 3958 expressed genes with P-value < 0.01 were screened out for subsequent TF and miRNA analyses.

Enrichment and annotation analysis results

The up- and down-regulated DEGs were inputted into the DAVID database to identify significantly enriched GO functional terms and KEGG pathways, respectively. Accordingly, the up-regulated DEGs were significantly enriched in biological processes related to pigmentation, such as the pigmentation during development (P-value = 3.63E-08), pigmentation [P-value = 5.78E-06; involving DCT (dopachrome tautomerase), SOX10 (SRY (sex determining region Y)-box 10), TYRP1 (tyrosinase-related protein 1), TYR (tyrosinase), MLPH (melanophilin), KIT (v-kit Hardy-Zuckerman 4 feline sarcoma viral oncogene homolog) and GPR143 (G protein-coupled receptor 143)], melanin biosynthetic process (P-value = 1.59E-03) and melanin metabolic process (P-value = 2.03E-03, involving DCT, TYRP1 and TYR) (Table 1).

Table 1

The significantly enriched GO_BP terms (top 1) and KEGG pathways of up-regulated DEGs

Term	Count	Gene lists	P-Value	Adjusted P-value
GO_BP: enrichment Score: 3.4367879252049858
GO_BP: 0048066 pigmentation during development	7	DCT, SOX10, TYRP1, TYR, MLPH, KIT, GPR143	3.63E-08	5.65E-05
GO_BP: 0043473 pigmentation	7	DCT, SOX10, TYRP1, TYR, MLPH, KIT, GPR143	5.78E-06	0.009009
GO_BP: 0042438 melanin biosynthetic process	3	DCT, TYRP1, TYR	1.59E-03	2.449285
GO_BP: 0006582 melanin metabolic process	3	DCT, TYRP1, TYR	2.03E-03	3.123208
GO_BP: 0046148 pigment biosynthetic process	4	DCT, TYRP1, TYR, GPR143	3.30E-03	5.019067
GO_BP: 0042440 pigment metabolic process	4	DCT, TYRP1, TYR, GPR143	4.96E-03	7.448208
GO_BP: 0030318 melanocyte differentiation	3	SOX10, TYRP1, MLPH	5.02E-03	7.534839
GO_BP: 0050931 pigment cell differentiation	3	SOX10, TYRP1, MLPH	5.76E-03	8.604376
KEGG
hsa04080: Neuroactive ligand-receptor interaction	6	SSTR4, GPR156, KISS1R, OPRL1, NPBWR1, GRIN1	1.58E-02	13.80265
hsa04916: Melanogenesis	4	DCT, TYRP1, TYR, KIT	1.99E-02	17.11413
hsa00350: Tyrosine metabolism	3	DCT, TYRP1, TYR	2.73E-02	22.81886

Note: GO gene ontology, BP biological process, KEGG Kyoto encyclopedia of genes and genomes, DEGs differentially expressed genes. Count, the number of DEGs enriched in the corresponding term. The gene symbols were listed in accordance with the Gene database at NCBI (National Center for Biotechnology Information, http://www.ncbi.nlm.nih.gov/gene/?term)

On the other hand, the down-regulated DEGs were dramatically related to the haemopoiesis progress and immune system, such as the haemopoiesis (P-value = 1.38E-04), haemopoietic or lymphoid organ development (P-value = 2.34E-04), immune system development [P-value = 3.22E-04; involving GPR183 (G protein-coupled receptor 143), INHBA (inhibin, beta A), PTPRC (protein tyrosine phosphatase, receptor type, C), PLEK (pleckstrin), CD8A (CD8a molecule) and IKZF1 (IKAROS family zinc finger 1)] (Table 2).

Table 2

The significantly enriched GO_BP terms (top 1) and KEGG pathways of down-regulated DEGs

Term	Count	Gene lists	P-Value	Adjusted P-value
GO_BP: enrichment Score: 2.505066250961214
GO_BP: 0030097 hemopoiesis	7	GPR183, INHBA, PTPRC, PLEK, CD8A, IKZF1, FLT3	1.38E-04	0.218579
GO_BP: 0048534 hemopoietic or lymphoid organ development	7	GPR183, INHBA, PTPRC, PLEK, CD8A, IKZF1, FLT3	2.34E-04	0.369612
GO_BP: 0050865 regulation of cell activation	6	INHBA, PTPRC, PLEK, IKZF1, FLT3, IL1B	3.05E-04	0.482267
GO_BP: 0002520 immune system development	7	GPR183, INHBA, PTPRC, PLEK, CD8A, IKZF1, FLT3	3.22E-04	0.509185
GO_BP: 0030098 lymphocyte differentiation	5	GPR183, PTPRC, CD8A, IKZF1, FLT3	4.05E-04	0.638643
GO_BP: 0002521 leukocyte differentiation	5	GPR183, PTPRC, CD8A, IKZF1, FLT3	1.00E-03	1.572197
GO_BP: 0045619 regulation of lymphocyte differentiation	4	INHBA, PTPRC, IKZF1, FLT3	1.27E-03	1.98708
GO_BP: 0045577 regulation of B cell differentiation	3	INHBA, PTPRC, FLT3	1.32E-03	2.064808
GO_BP: 0030217 T cell differentiation	4	PTPRC, CD8A, IKZF1, FLT3	1.39E-03	2.173407
GO_BP: 0051249 regulation of lymphocyte activation	5	INHBA, PTPRC, IKZF1, FLT3, IL1B	1.57E-03	2.458734
KEGG pathways
hsa04640: Hematopoietic cell lineage	7	CD8A, FLT3, CD1C, IL1B, CSF2RA, CD1E, IL1A	1.18E-06	0.001103
hsa04060: Cytokine-cytokine receptor interaction	8	INHBA, CCL22, FLT1, IL18RAP, FLT3, IL1B, CSF2RA, IL1A	7.78E-05	0.072785

Note: GO gene ontology, BP biological process, KEGG Kyoto encyclopedia of genes and genomes, DEGs differentially expressed genes. Count, the number of DEGs enriched in the corresponding term. The gene symbols were listed in accordance with the Gene database at NCBI (National Center for Biotechnology Information, http://www.ncbi.nlm.nih.gov/gene/?term)

Additionally, functional annotation analysis revealed that the down-regulated DNAJB4 (DnaJ (Hsp40) homolog, subfamily B, member 4) and SLIT2 (slit homologue 2) might function as tumor-suppressors, and the up-regulated KIT (v-kit Hardy-Zuckerman 4 feline sarcoma viral oncogene homolog) may function as an oncogene.

Protein-protein interaction network

By inputting the DEGs into the STRING database, a total of 152 interaction pairs was screened out and then constructed into a PPI network (Fig. 1). In the network, the top 3 hub nodes with highest connectivity degree were PTPRC (degree = 13), TYR (degree = 13), and PLEK (degree = 12).

Fig. 1

Protein-protein interaction network of differentially expressed genes. Circular nodes represent differentially expressed genes and edges represent interactions. Red represents up-regulation while green represents down-regulation. The color depth represents the statistical significance level

Regulatory networks

Furthermore, a number of TFs and miRNAs were screened from the DEGs. From the Fig. 2, significantly up-regulated KIT is a potential target of miR-206 and miR-496, and DCT may be targeted by miR-184. MiR-337-5p (P-value = 7.62E-04) was most significantly enriched for the up-regulated DEGs while miR-21 (P-value = 3.36E-04) and miR-16 (P-value =8.90E-04) were identified as having the most significance for down-regulated DEGs (Table 3). Moreover, based on Fig. 3, the up-regulated PAX5 (paired box 5) and HNF4A (hepatocyte nuclear factor 4, alpha) targeted MLPH and GPR143, respectively, while the down-regulated BATF (basic leucine zipper transcription factor, ATF-like), SPI1 (Spi-1 proto-oncogene) and EP300 (E1A binding protein p300) jointly targeted GPR183, PTPRC and PLEK. The activity differences between UVB-irradiated and non-irradiated skin biopsies were identified as 4.78E-05, 5.24E-02, −1.06E-01, −1.05E-01 and −8.42E-02 for PAX5, HNF4A, BATF, SPI1 and EP300, respectively (Table 4).

Fig. 2

Regulatory network of microRNA-target genes. Rhombuses represent microRNAs and circles represent differentially expressed genes. Red represents up-regulation while green represents down-regulation. –—, interactions between DEGs; −—|, regulatory effects of microRNAs on DEGs. DEGs, differentially expressed genes

Table 3

The miRNAs targeting the differentially expressed genes (DEGs)

MicroRNA	Count	P-value
Targeting up-regulated DEGs
hsa-miR-337-5p	4	7.62E-04
hsa-miR-326	5	1.50E-03
hsa-miR-496	4	1.66E-03
hsa-miR-515-3p	3	2.53E-03
hsa-miR-542-5p	3	5.91E-03
hsa-miR-510	2	7.75E-03
hsa-miR-142-3p	9	8.97E-03
hsa-miR-142-5p	7	1.07E-02
hsa-miR-337-3p	4	1.18E-02
hsa-miR-133b	10	1.69E-02
hsa-miR-375	9	1.87E-02
hsa-miR-206	11	1.90E-02
hsa-miR-299-3p	3	2.71E-02
hsa-miR-455-5p	3	2.71E-02
hsa-miR-429	10	3.02E-02
hsa-miR-184	5	3.20E-02
hsa-miR-339-5p	3	3.81E-02
hsa-miR-455-3p	3	3.81E-02
Targeting down-regulated DEGs
hsa-miR-21	12	3.36E-04
hsa-miR-16	20	8.90E-04
hsa-miR-15a	18	1.55E-03
hsa-miR-296-3p	9	5.86E-03
hsa-miR-106a	3	1.07E-02
hsa-miR-135a	2	2.09E-02
hsa-miR-129-5p	5	4.00E-02
hsa-miR-342-3p	5	4.00E-02
hsa-miR-342-5p	5	4.00E-02
hsa-miR-296-5p	7	4.71E-02

Fig. 3

Regulatory network of transcription factors-target genes. Triangles represent transcription factors and circles represent differentially expressed genes. Red represents up-regulation while green represents down-regulation. –—, interactions between differentially expressed genes; →, regulatory effects of transcription factors on differentially expressed genes

Table 4

The activity differences of transcription factors

Transcription factors	Activity difference
Up-regulated
NFE2	3.68E-01
ZBTB33	3.28E-01
NR3C1	2.82E-01
HDAC8	2.34E-01
SMARCC2	1.34E-01
CCNT2	1.10E-01
PPARGC1A	5.60E-02
HNF4A	5.24E-02
GATA2	2.12E-02
PAX5	4.78E-05
Down-regulated
BRF1	−3.95E-01
CTCF	−2.34E-01
YY1	−1.68E-01
ESRRA	−1.48E-01
BATF	−1.06E-01
SPI1	−1.05E-01
ELF1	−9.63E-02
EP300	−8.42E-02
E2F6	−6.76E-02
REST	−1.88E-02
PBX3	−1.49E-02

Note: Activity difference = Activity in UVB-irradiated skin biopsies – Activity in non-irradiated skin biopsies

Integrated network and screened modules

By merging the PPI network and two regulatory networks, an integrated network among TFs, miRNAs and DEGs was obtained (Fig. 4), which was then mined to obtain two modules with highest score (Fig. 5). In the first module, there were 9 down-regulated DEGs including one TF, EP300, which is significantly associated with the immune system, such as the immune response (involving GPR183 and PTPRC), haemopoiesis, haemopoietic or lymphoid organ development, immune system development and cell activation (involving GPR183, PTPRC and PLEK) (Table 3). In contrast, there were 7 up-regulated DEGs in the second module, which are dramatically related to pigmentation, including the pigmentation during development, pigmentation, pigment biosynthetic process and pigment metabolic process (involving DCT, TYRP1, TYR and GPR143) (Table 5).

Fig. 4

Integrated network among transcription factors, microRNAs and differentially expressed genes. Rhombuses represent microRNAs, triangles represent transcription factors and circles represent differentially expressed genes. Red represents up-regulation while green represents down-regulation. –—, interactions between differentially expressed genes; −—|, regulatory effects of microRNAs on differentially expressed genes; →, regulatory effects of transcription factors on differentially expressed genes

Fig. 5

Two densely connected modules screened from the integrated network. a down-regulated module; b up-regulated module. Triangles represent transcription factors and circles represent differentially expressed genes. Red represents up-regulation while green represents down-regulation. –—, interactions between differentially expressed genes; →, regulatory effects of transcription factors on differentially expressed genes

Table 5

The top 5 significantly enriched GO_BP terms of screened modules

Term	Count	Gene lists	P-Value	Adjusted P-value
Module 1 Score: 7.75
GO_BP: 0006955 ~ immune response	4	GPR183, PTPRC, NCF4, GPR65	2.35E-03	3.251619
GO_BP: 0030097 ~ hemopoiesis	3	GPR183, PTPRC, PLEK	4.34E-03	5.921337
GO_BP: 0048534 ~ hemopoietic or lymphoid organ development	3	GPR183, PTPRC, PLEK	5.25E-03	7.113787
GO_BP: 0002520 ~ immune system development	3	GPR183, PTPRC, PLEK	5.89E-03	7.959213
GO_BP: 0001775 ~ cell activation	3	GPR183, PTPRC, PLEK	6.36E-03	8.562876
Module 2 Score: 6.333
GO_BP: 0048066 ~ pigmentation during development	5	DCT, SOX10, TYRP1, TYR, GPR143	1.07E-11	1.10E-08
GO_BP: 0043473 ~ pigmentation	5	DCT, SOX10, TYRP1, TYR, GPR143	3.26E-10	3.35E-07
GO_BP: 0046148 ~ pigment biosynthetic process	4	DCT, TYRP1, TYR, GPR143	8.84E-08	9.07E-05
GO_BP: 0042440 ~ pigment metabolic process	4	DCT, TYRP1, TYR, GPR143	1.37E-07	1.41E-04
GO_BP: 0042438 ~ melanin biosynthetic process	3	DCT, TYRP1, TYR	1.84E-06	0.001882

Note: GO gene ontology, BP biological process. Count, the number of differentially expressed genes enriched in the corresponding term. The gene symbols were listed in accordance with the Gene database at NCBI (National Center for Biotechnology Information, http://www.ncbi.nlm.nih.gov/gene/?term)

Discussions

Using expression profiling of UV radiation-induced melanocytes, Yang et al. revealed a set of p53 target genes in response to UV radiation through the Ingenuity® Pathway Analysis program [27]. In 2015, López et al. conducted a time-series analysis to compare the transcriptional profile of dark and light melanocyte lines under basal conditions and after UVB radiation, finding that ribosomal proteins can interact with the p53 signaling pathway in response to UVB [28]. This study identified a total of 151 up-regulated and 64 down-regulated genes in skin in response to UVB irradiation, which probably remained increased after repetitive exposure for approximately 2–3 weeks [11]. The up-regulated DEGs were significantly related to pigmentation, whereas the down-regulated DEGs were dramatically associated with immunity. The functional enrichment results were observed and were similar in the screened modules. As mentioned above, the up-regulated DEGs were significantly enriched in biological processes related to pigmentation. The obtained results are consistent with a previous study [11], implying the reliability of our research. However, the p53 pathway was not enriched in this study, indicating that in vitro and in situ studies have discrepancies and that an in situ study might be more credible. Importantly, this study further identified a number of TFs and miRNAs potentially targeting DEGs.

Several DEGs were enriched in pigmentation, such as DCT, SOX10, TYRP1, TYR, MLPH, KIT and GPR143, whose relationships with pigmentation have been revealed extensively in a previous study [11]. However, this study further revealed close interactions among DCT, SOX10, TYRP1, TYR and GPR143 in a densely connected module, suggesting a collaborative stimulating effect on the pigmentation process. Moreover, TYR had many interactions with other DEGs and was a hub node in the PPI network, implying predominant roles in the biological processes of pigmentation. Additionally, the up-regulated KIT, which was previously identified as a oncogenic drivers in sun-impaired skin, was annotated as a potential oncogene in this study [29].

Notably, some TFs and miRNAs were screened from the up-regulated DEGs. KIT is a potential target of miR-206 and miR-496, and DCT may be targeted by miR-184. As previously reported, miR-206 acts as a tumor suppressor in melanoma via inhibiting cancer cell growth and migration [30], miR-496 is involved in silencing gene expression in breast cancer, and miR-184 antagonizes miR-205 to protect against stratified squamous epithelia [31]. In consideration of the involvement of KIT and DCT in pigmentation and carcinogenesis, miRNAs might also be implicated in pigmentation in response to UVB irradiation. Moreover, miR-337, which can have a tumor suppressive effect on pancreatic ductal adenocarcinoma [32], was most significantly enriched, implying its involvement in the cellular response to UVB exposure. On the other hand, PAX5 and HNF4A were identified as two TFs targeting MLPH and GPR143, respectively. MLPH has been demonstrated as a putative pigmentation genes in East Asia [33], and GPR143 is targeted to melanosomes as a G protein-coupled receptor in pigment cells [34]. Thus, PAX5 and HNF4A might be involved in the regulation of pigmentation under UVB exposure. The up-regulation of DEGs may be associated with the screened the miRNAs and TFs under UVB irradiation.

On the other hand, the down-regulated DEGs were related to the immune system, suggesting an inactive immune regulation, e.g., GPR183, INHBA, PTPRC, PLEK, CD8A and IKZF1. UVB irradiation is responsible for inducing systemic immune suppression, which is an important contributor to skin cancer [35]. The relevance of those DEGs to the immune system or cancer has been revealed. Constitutively activated PTPRC modulates neutrophil recruitment and integrin activation during inflammation [36]. GPR183, which is primarily expressed in lymphoid cells, is important for the adaptive immune response [37]; as well as PLEK and CD8A are also important in this response [38, 39]. In addition, IKZF1 encodes a TF regulating lymphocyte differentiation [40], and INHBA, which belongs to the TGF-β superfamily, is associated with several cancers [41]. In the screened module, we further identified close interactions among GPR183, PTPRC and PLEK and observed their enrichment in the immune system. Notably, PTPRC and PLEK were hub nodes in the PPI network. Therefore, these down-regulated DEGs might jointly contribute to the suppressive immune response to UVB exposure, which might be a cause of skin cancer development. Additionally, the down-regulated DNAJB4, which is implicated in human gastric carcinomas [42], and SLIT2, which is involved in pancreatic cancer [43], are also annotated as tumor suppressors, suggesting their potential roles in UVB-induced skin cancers.

Additionally, some miRNAs were identified for the down-regulated DEGs, among which miR-21 and miR-16 were the most significantly enriched. The expression of miR-21 could be decreased by UVB irradiation in NIH3T3 cells [44]. MiR-16 targets the thyroid hormone receptor and nuclear factor-kappaB, which are vital for initiating epithelial immune responses [45]. Thus, miR-21 and miR-16 may play a role in the UVB-induced immunosuppressive response via down-regulating immune-related genes. In addition, GPR183, PTPRC and PLEK are all targeted by TFs BATF, SPI1 and EP300, which have been demonstrated to be involved in the immune system [46-48]. Notably, we also observed interactions of EP300 with other down-regulated DEGs in the screened module. The p53/p300 complex can be recruited to the promoters of p53 target genes that are related to cell growth and apoptosis under UVB irradiation [49]. Therefore, the TFs and miRNAs may be involved in the cellular immune response to UVB exposure.

Conclusion

In conclusion, this study identified 151 up-regulated and 64 down-regulated genes implicated in the pigmentation and immune system response to UVB exposure. UVB-induced pigmentation has been demonstrated in previous studies, but the down-regulation of respectively by UVB exposure involved in the immune system have not been revealed. Furthermore, this study screened some TFs and miRNAs for the DEGs, which might play important roles in UVB irradiation-induced pigmentation and the immunosuppressive response. However, although we identified many DEGs, TFs and miRNAs with potential critical roles in UVB irradiation-induced cellular response, the obtained results relied on bioinformatics methods and need further validation.