The microarray and RNA-sequencing (RNA-Seq) datasets were downloaded from the Gene Expression Omnibus¹ and Sequence Read Archive² databases to obtain gene expression profiles using the keywords “UVB” and “UVB radiation.” Basic inclusion criteria were (1) gene expression profiles of human-derived skin samples (2). Experiments that exposed samples to UVB without other substances or drugs (3). Acute exposure experiments, which refer to contact with a substance that occurs once or for only a short time (14). Finally, (4) studies that provided processed RNA-Seq data where read count data must be accessible to prevent noise from the analysis results. The final selected datasets consisted of two microarray experiments and six RNA-Seq experiments, encompassing 39 control samples and 41 UVB-exposed samples. These datasets and their detailed information are summarized in Supplementary Table S1.

The raw expression data of the Agilent microarray (GSE45493) and Affymetrix (GSE41078) microarray were normalized using the “normalizedBetweenArrays” function in the limma R package (15).

In the RNA-Seq data, raw data (FASTQ files) were downloaded from the SRA database, and low-quality reads were trimmed using Trimmomatic software (v0.36) (16). The trimmed reads were aligned to the human reference genome (GRCh38) using HISAT2 (v2.2.1) (17). SAMtools (v1.10) was used to convert the SAM file to a BAM file and sort the sequences (18). The read count data were generated by StringTie (v2.2.0) or downloaded from the GEO database (19). Low-count genes were filtered out, and normalization was performed using edgeR packages (v3.38.1) in R software (20). The expression level of RNA-Seq results was analyzed using counts per million (CPM) values.

To integrate microarray and RNA-Seq data, the CPM values of RNA-Seq data were transformed (21, 22) using Eq. (1) to eliminate the singularity at CPM = 0.

After transformation, each dataset was annotated with an Ensembl ID and expression data were integrated based on common Ensembl gene IDs. Batch effects (non-biological differences) were adjusted by applying the “sva” function (v3.44.0) in R software (23).

After batch-effect correction, the principal component analysis (PCA) was conducted to examine the distributions derived from different groups. Differentially expressed genes (DEGS) between the control and UVB-exposed groups were identified using the limma R package (v3.52.2). The threshold criteria were |fold change| ≥1.5, and false discovery rate adjusted p value < 0.05. The expression patterns of the DEGs were visualized using a hierarchical clustering heatmap.

Pathway Studio (v12.5.0.2), a commercial text mining-based analytical tool, was used to explore and visualize the various types of functional interactions between biological entities, including genes, cell processes, and diseases (24). Pathway Studio provides insight into the biological relationships between the entered genes, as well as an investigation of protein-functional class, protein-disease, or protein-cell process. In addition, to uncover the key entities that contribute significantly to the network, betweenness centrality and degree values were calculated using NetworkAnalyzer in Cytoscape (v3.8.2) (25). To construct a network based on more genetically significant genes, threshold criteria were narrowed (|fold change| ≥2). Finally, the top 45 genes were selected with a degree ≥3 based on the gene–gene interactions derived from the Pathway Studio database information. In addition, the signaling network was expanded to explore major cell processes, cutaneous diseases, and functional classes associated with selected genes. Twenty-nine cell processes, 21 skin diseases, and 5 functional classes were selected based on the criteria that the number of references should be at least 3 for all relations, and the degree value should be at least 3 for skin disease, 5 for cell process, and 10 for functional class.

To explore the potential function of the identified significant DEGs, gene ontology term analysis was performed using the clusterProfiler R package (v4.4.4) (26). The gene ontology or GO terms, refer to the following domains: biological process, cellular component, and molecular function. Enrichment analysis was performed with the threshold of adjusted p value < 0.05.

A survival analysis was performed using Survival Genie (https://bbisr.shinyapps.winship.emory.edu/SurvivalGenie/), a web platform for survival analysis related to adult and pediatric cancers (27). Genes were uploaded to the Survival Genie web platform, and SKCM dataset in The Cancer Genome Atlas (TCGA) dataset was selected to perform the survival analysis. Metastatic tumor type and median-based partitioning method were selected. Kaplan–Meier survival curves were generated by the Survival Genie platform, and the survival rates of the stratified high (red) and low (blue) groups of patients were shown according to log-rank tests.

Normal Human Epidermal Keratinocytes (NHEKs) and Normal Human Dermal Fibroblasts (NHDFs) were purchased from Lonza (Basel, Switzerland). NHEKs were cultured in keratinocyte growth medium (KBM-Gold) with SingleQuots™ supplement (Lonza), including bovine pituitary extract, epinephrine, gentamicin/amphotericin B, human epidermal growth factor, hydrocortisone, insulin, and transferrin at 37°C in a 5% CO₂ incubator. NHDFs were cultured in Dulbecco’s modified Eagle medium (DMEM) containing 4.5 g L⁻¹ glucose, L-glutamine, 10% fetal bovine serum, and 1% penicillin–streptomycin mixture (Gibco, United States) at 37°C in a 5% CO₂ incubator. NHEKs and NHDFs were irradiated with UVB (312 nm) at various intensities (11–20 mJ cm⁻²) in phosphate-buffered saline using a Bio-Sun UV irradiation system (Vilber Lourmat). After UVB exposure, the NHEKs were incubated in KBM-Gold and the NHDFs were incubated in DMEM, both for 24 h.

Cell viability was measured to select an appropriate UVB intensity. NHEKs were irradiated at 13 mJ cm⁻² based on the results of the cell viability test in our previous study (28). For NHDFs, the 3-(4,5-dimethylthiazol-2-yl)-2,5 diphenyltetrazolium bromide (MTT) assay was performed to confirm cell viability. NHDFs were seeded in 96-well cell culture plates and stabilized overnight at 37°C. NHDFs were shaded with aluminum foil and irradiated with 0, 11, 13, 15, and 19 mJ cm⁻² UVB, respectively. After 24 h at the end of the exposure, the MTT assay was performed according to the manufacturer’s protocols.

RNA from NHEKs and NHDFs was extracted and the quality control of RNA samples was performed as described previously (28). cDNA was synthesized and quantitative reverse transcription polymerase chain reaction (qRT-PCR) was performed as described previously (28). The primer sequences of four genes, GAPDH, IL1B, CCL2, and LIF, and the thermal cycling conditions are shown in Supplementary Table S2.

Statistical analyses were performed using R software (v4.2.2). The differences between the control and exposure groups were analyzed using Student’s t-test with p-value <0.05. GO analyses were performed using R software, and adjusted p-value <0.05 were considered statistically significant. Statistical analysis of survival analysis was performed in Survival Genie software, and a log-rank p-value <0.05 was considered significant.

To identify consistently regulated genes upon UVB exposure of skin, eight gene expression profile datasets were collected from the GEO and SRA databases considering the inclusion criteria. Specifically, six different RNA-Seq datasets were used, including one previously published dataset (28) and two microarrays. Furthermore, various skin cell types, such as skin biopsies, fibroblasts, and keratinocytes, were used to prevent biases toward single cell types. To overcome differences between the microarray and RNA-Seq platforms, the RNA-Seq data were slightly transformed (as detailed in the “Materials and methods”). As the merged expression dataset originated from various cell types and platforms, it was important to focus solely on the effect of UVB exposure. Therefore, the batch effects were adjusted for each cell type to focus on UVB exposure, as shown in Figure 1A.

The results of PCA indicated that UVB-exposed samples were mostly distinct from control samples with some exceptions, and most of the variances due to tissue origin or the platform of each dataset were not identified (Figure 1B). After the integration of eight studies, 681 robust genes that were differentially expressed were extracted (upregulated: 450, downregulated: 231; Figure 1C). There was an almost complete consistent distinction in the gene expression profiles between the control and UVB exposure groups (Figure 1D).

GO enrichment analysis was performed to investigate the potential functions of the obtained DEGs. Figure 1E and Supplementary Table S3 show which biological functions were mainly predicted by upregulated genes. In biological processes, the activated GO terms were confirmed to be related to the immune response, such as response to lipopolysaccharide, cytokine-mediated signaling pathway, and response to a molecule of bacterial origin. In cellular components, the main GO terms were related to intracellular organelles, including tertiary granule, ficolin-1-rich granule, vesicle lumen, secretory granule lumen, and cytoplasmic vesicle lumen. In molecular functions, cytokine responses, including cytokine activity and cytokine receptor binding, were the predicted enriched GO terms, and reactions to temperatures and external stimuli, including heat shock protein binding, were identified.

At the same time, we identified biological functions that were mainly predicted by downregulated genes (Figure 1F; Supplementary Table S4). In biological processes, the suppressed GO terms were confirmed to be related to xenobiotic response, including xenobiotic metabolic process, cellular response to xenobiotic stimulus, and response to xenobiotic stimulus. In cellular component, there was only one related GO term; the downregulated genes were mainly involved in the complex of collagen trimers. In molecular functions, the downregulated genes were mostly associated with two categories: (1) oxidoreductase activity, acting on paired donors, with the incorporation or reduction of molecular oxygen, and (2) glutathione transferase activity.

To explore the dermatologic signaling network and hub genes associated with UVB exposure, we used Pathway Studio software to perform network analysis of DEGs in response to UVB exposure. Based on the text-mining algorithm of the software, we expanded the signaling network to predict the biological functions that may be altered in UVB-exposed skin (Figure 2A). The genes significantly associated with inflammatory signaling were identified from functional classes, including NF-κB family, ERK1/2, mitogen-activated protein kinase, JNK, and Jun/Fos. In terms of cell processes, the network among the genes was closely related to altered carcinogenic processes, including tumor growth, cell invasion, and cancer cell growth. Inflammatory cell processes were also predicted, including inflammatory response and immune response. Association with cancer-related diseases was shown, including skin cancer, basal cell carcinoma, metastatic melanoma, and nonmelanoma skin cancer.

We analyzed the biological hub network from Figure 2A to construct potential hub networks that might elucidate the biological relationships between cutaneous effects and UVB exposure (Figure 2B). Seven hub genes (IL1B, IL6, MMP9, CCL2, CXCL8, PTGS2, and LIF) and major entities were screened considering the betweenness centrality and degree value (Supplementary Table S5; Supplementary Figure S1). NF-κB, Jun/Fos, and mitogen-activated protein kinase were predicted to be the major functional classes. Tumor growth, inflammatory response, immune response, wound healing, and apoptosis were predicted to be the major cellular processes. In major diseases, psoriasis, dermatitis, vitiligo, skin cancer, metastatic melanoma, and nonmelanoma skin cancer were selected.

We performed a survival analysis to examine the correlation between the relative expressions of the seven hub genes and patients with metastatic SKCM. All patients with metastatic SKCM were divided into low- and high-expression groups according to the gene expression levels (Table 1). For the set of 7 hub genes, 179 low- and 178 high-expression groups were separated among 357 patients with SKCM. As shown in Figure 3A, the set of seven hub genes was considered to be significantly associated with the prognosis (p value < 0.05) of patients with SKCM. We found that IL1B, CCL2, and LIF among the seven hub genes were individually associated with the overall survival of patients with SKCM (p value < 0.05; Table 1; Figure 3A).

To confirm the expression level of the three hub genes (IL1B, CCL2, and LIF) in keratinocytes and fibroblasts upon UVB exposure, we used NHDFs and NHEKs. The UVB exposure resulted in dose–response curves for both cell lines. UVB intensity showed >80% cell viability at an irradiation dose of 13 mJ cm⁻² in NHEKs and 11 mJ cm⁻² in NHDFs (Supplementary Figure S2). qRT-PCR was performed to confirm the RNA expression levels of hub genes under UVB exposure conditions. Among them, all three hub genes showed significant upregulation in UVB-irradiated NHDFs and NHEKs (Figures 3B,C).