Identification and validation of four photodynamic therapy related genes inhibiting MAPK and inducing cell cycle alteration in squamous cell carcinoma.

Introduction: Cutaneous squamous cell carcinoma (cSCC) is the second most common skin cancer, and photodynamic therapy (PDT) is a promising modality against cSCC. This study investigated the impact of PDT on the MAPK pathway and cell cycle alternation of cSCC as well as the related molecular mechanisms. Method: Expressing mRNA profile data sets GSE98767, GSE45216, and GSE84758 were acquired from the GEO database. The functions of differently expressed genes (DEGs) were enriched by Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. Least absolute shrinkage and selection operator (Lasso) analysis were used to establish a diagnosis model based on GSE98767. A correlation analysis and a protein-protein interaction (PPI) network were used to evaluate the relationship between cSCC-PDT-related genes and the MAPK pathway. Single-sample gene set enrichment analysis (ssGSEA) was performed on GSE98767 to estimate MAPK activation and cell cycle activity. Finally, the effect of MAPK activation on the cell cycle was explored in vitro. Result: Four cSCC-PDT-related genes, DUSP6, EFNB2, DNAJB1, and CCNL1, were identified as diagnostic markers of cSCC, which were upregulated in cSCC or LC50 PDT-protocol treatment and negatively correlated with the MAPK promoter. Despite having a smaller MAPK activation score, cSCC showed higher cell cycle activity. The PDT treatment suppressed the G1 to G2/M phase in JNK overexpressed A431 cells.
Conclusion: CCNL1, DNAJB1, DUSP6, and EFNB2 were identified as potential PDT target genes in cSCC treatment, whose potential therapeutic mechanism was inhibiting the MAPK pathway and inducing cell cycle alternation.

Introduction

Cutaneous squamous cell carcinoma (cSCC), a malignant tumor originating in epithelial cells, is the second most common skin cancer, with an increasing incidence in elderly individuals (1). Surgical excision is the first-line treatment for common primary cSCC (2), but the obvious scars resulting from extended resection, especially those on eyelids, ear lesions, or lips, lead to poor cosmetic outcomes. Moreover, patients on immune suppressants or anticoagulants and elderly individuals with comorbid conditions are ineligible for surgery. Nonsurgical approaches such as radiotherapy, chemotherapy, and cryotherapy are applied to the aforesaid patients. Unfortunately, these treatments are associated with severe side effects, chemoresistance, or a high recurrence rate (3). Thus, there is an urgent need to develop a novel therapy with a minimally invasive nature, good anti-tumor effect, and limited side effects.

Photodynamic therapy (PDT) is a burgeoning, safe, effective, and affordable remedy that provides a great cosmetic outcome for cSCC. PDT requires three components: a photosensitizer, light, and molecular oxygen. In the treatment process, the photosensitizer is first coated or injected into the tumor tissue, and then the nidus is irradiated by the light at the appropriate wavelength and energy. The photosensitizer in tumor tissues can produce a photochemical reaction to generate reactive oxygen species (ROS), which will react with bioactive components in tumor cells to kill them (4). Additionally, ROS can be rapidly metabolized in normal tissues around tumors, which signifies that PDT is highly selective for treating cSCC.

Previous mechanism-based studies have revealed that PDT induces apoptosis or necrosis of cSCC cells primarily through the direct cytotoxicity of ROS, complete occlusion of microvasculature that nourishes cancer cells, and stimulation of anti-tumor immunity (5). Furthermore, increasing evidence has shown that the occurrence and progression of cSCC may be related to abnormal regulation of mitogen-activated protein kinase (MAPK) signaling. The mutation of the tumor suppressor genes TP53, CDKN2A, NOTCH, and oncogene Ras is common in cSCC (6–8). The mutation of Ras leads to the declining activity of GTPase-activating proteins (GAPs), suppressing the transformation from active guanosine triphosphate (GTP)-bound states to inactive guanosine diphosphate (GDP)-bound states (9). Redundant GTP causes Ras activation and upregulation of downstream PI3K/AKT/mTOR and MAPK intracellular signaling (10) and eventually leads to excessive proliferation of cells and tumorigenesis. The best-known MAPKs include c-Jun amino-terminal kinases 1 to 3 (JNK1 to 3), extracellular signal-regulated kinases 1 and 2 (ERK1/2), ERK5 and p38 (α, β, γ, and δ) families (11). ERK is related to cell growth and differentiation, while JNK is related to inflammation and apoptosis. Previous studies (12) have proved that the suppression of JNK-dependent apoptosis, cooperating with paradoxical ERK activation, induces cSCC. To sum up, the progression of cSCC is closely associated with MAPK signaling (13). Moreover, some studies (14, 15) indicate that MAPK activation stimulates cell cycle transitions, and numerous therapies have exerted their anti-tumor effects by blocking the MAPK signaling pathway and thus arresting the cell cycle transformation. A recent study (16) reported that modified 5-aminolevulinic acid photodynamic therapy could suppress the ERK pathway, activate the JNK pathway, and inhibit cell proliferation and Cyclin-D1, which performed a regulatory function during the transformation from G1 to G2/M phase. Accordingly, we can infer that PDT probably affects cell cycle alterations in cSCC through the MAPK pathway.

Preceding bioinformatics analysis suggested that the combination of PDT and chemotherapy could upregulate ROS generation, blocking the cell cycle and proliferation (17). Besides, it seems that PDT downregulates the pathways related to tumor progression, metastasis, and invasion. Up to now, there are still rare reports about the comprehensive mechanism of PDT for cSCC treatment. We were interested in exploring whether inhibition of cSCC proliferation by PDT was related to MAPK signaling and its regulation of the cell cycle. In this study, we found a potential molecular mechanism of PDT for treating cSCC was that PDT inhibited MAPK activation, which was probably because the downregulation of JNK was greater than the upregulation of ERK. We preliminarily confirmed that the inhibition of MAPK suppressed the transformation from G1 to G2/M phase in cSCC using bioinformatics analysis and flow cytometry. Additionally, CCNL1, DNAJB1, DUSP6, and EFNB2 were upregulated in cSCC or LC50 PDT-protocol treatment based on our results. Thus, we speculated that increased CCNL1, DNAJB1, DUSP6, and EFNB2 inhibited the MAPK pathway in cSCC, and they were potential PDT target genes in cSCC treatment.

Method

Data sources

Expressing mRNA profiles and associated clinical data of cSCC were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo/). Dataset GSE98767, which was performed on GPL10558 and consisted of 15 cSCC patients and three normal subjects with three repetitions (54 samples), was used as the training set to screen different expression genes (DEGs) (18). Dataset GSE45212 performed on GPL570 was used as the validating set, which included 30 cSCC samples and 10 actinic keratosis (AK) samples (40 samples) (19). Dataset GSE84758 was performed on GPL10558 and contained 48 A431 cell samples, among which three control samples, three DT samples, three LC50 PDT-protocol treated samples, and three LC90 PDT-protocol treated samples (12 samples) were used to filter DEGs after PDT treatment (20).

DEG analysis for normal and tumor groups

The “limma” R package was used to screen DEGs after log2(counts+1) standardizing. The P-value was adjusted by the Bayes test. The adjusted P-value<0.05 and |Log2 (Fold Change) | >1 were set as the cutoff, and the DEGs were applied to further analysis.

Functional analysis for DEGs

To assess the potential functions of DEGs, the “clusterProfiler” R package was applied to execute Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis. The P-value<0.05 was set as the cutoff criterion.

Least absolute shrinkage and selection operator regression establish the diagnostic model

The “glmnet” R package and least absolute shrinkage and selection operator regression (Lasso) analysis were used to build the cSCC diagnosis model, thus further filtering cSCC-PDT-related genes. The diagnostic score of each sample was calculated using the following formula: . The diagnostic score and expression level of selected genes with the maximum AUC value were selected as the boundary values, and the receiver operator characteristic (ROC) curve was drawn in both the training set and validating set.

Protein–protein interaction network

Protein–protein interaction (PPI) data were obtained from the STRING database (https://cn.string-db.org/, version 11.5) and processed in Cytoscape software.

Single sample gene set enrichment analysis

The Single sample gene set enrichment analysis (ssGSEA) was applied to explore the enriched scores of MAPK activation level, angiogenesis, cell cycle, differentiation, and proliferation in GSE98767 and GSE84758 using the R package “GSVA” (21). All used GO molecular functional (GOBP) gene sets were acquired from MSigDB (22) (https://www.gsea-msigdb.org/gsea/msigdb/index.jsp, version 7.5.1) and are annotated in .

Single-cell validation

The CancerSEA database (http://www.biocc.hrbmu.edu.cn/CancerSEA/home.jsp) (23) and EXP0063 were used to validate the correlation between the expression of selected genes and functional states in single-cell lines. EXP0063 (GSE103322) (24) was performed on GPL18573 containing RNA-seq data of 5,902 single cells from 18 head and neck squamous cell carcinoma (HNSCC) patients. The outcome of the tSNE plot and functional states were downloaded from the CancerSEA online tool.

Cell culture and cell cycle analysis

A431 cells (Thermo Fisher Scientific, USA) were cultured in the 96-well plates at a density of 5 × 103 cells/well and incubated for 24 h. The following compounds were used to treat A431 cells: the EGF group was treated with 100 ng/ml recombinant EGF (sc-4552, Santa Cruz Biotechnology, USA), ERK overexpression group was treated with 10 µM Ceramide C6 (sc-3527, Santa Cruz Biotechnology, USA), and the JNK overexpression group was treated with 500 nM anisomycin (sc-3524, Santa Cruz Biotechnology, USA). Based on our previous research, 400 nmol/L of EtNBSe (Department of Chemistry of Central South University, China) was used to pretreat all groups in the dark for 20 min. Then the cells were irradiated by LED light whose wavelength was 635 nm and power density of 50 mW/cm2 for 60 s (2.8 J/cm2) (25, 26). After corresponding treatments, cells were harvested and fixed with 75% ethanol overnight at 4°C. The fixed cells were washed three times with PBS and stained darkly with propidium iodide (PI) for 30 min. The stained cells were analyzed by a FACS Calibur flow cytometer (BD Biosciences, USA).

Statistical analysis

A Spearman analysis was performed to calculate the correlation coefficient between different genes. All statistical analyses were executed in the R language (version 4.1.3) and p<0.05 was considered statistically significant.

Result

The 38 cSCC-PDT related genes enriched in the MAPK pathway

The study design is illustrated in . The nine normal samples and 45 cSCC samples from GSE98676 were used as training sets for differential gene analysis. After Log2(counts+1) standardization, the signal density of each sample was generally at the same level for comparison ( ). Principal component analysis (PCA) showed that the two groups had different characteristics ( ). Following the analysis of the GSE98767 dataset using limma, a total of 2,293 DEGs were identified in the cSCC and normal samples, which were displayed in the volcano plot ( ). Dataset GSE84758 was processed in the same way as GSE98767, and 183 DEGs were identified. Then, 38 cSCC-PDT-related genes were obtained from the intersection of the above two DEG sets ( ). visually displays the expression level of the 38 genes in GSE98767, while lists information on the 38 genes in detail.

Figure 1

Identification of differentially expressed cSCC-PDT-related genes. (A) Box plot of 54 samples of GSE98767 after log2 standardization original counts. (B) Principal component analysis (PCA) analysis of tumor and a normal group of GSE98767. (C) Volcano plot of DEGs in cSCC and normal samples of GSE98767. The red square dots represented upregulated genes, while the green square dots represented downregulated genes (cutoff: adj. P-value<0.05 and |Log2 (Fold Change) |>1). (D) Venn diagram showing the 38 cSCC-PDT-related genes (the intersection of the cSCC and PDT). (E) Heatmap displaying the expression level of 38 cSCC-PDT-related genes in GSE 98767.

Table 1

Detailed data of 38 differential expressed cSCC-PDT-related genes.

Gene.Symbol	GSE98767			GSE84758
	logFC	P.Value	adj.P.Val	logFC	P.Value	adj.P.Val
HIST1H2BK	4.619604846	3.2167E-17	1.87778E-16	1.31911	6.12E-05	0.009976
ERRFI1	4.277016012	1.41405E-20	9.2757E-20	2.833846	1.29E-08	1.88E-05
DUSP6	4.032873036	8.16355E-30	7.60036E-29	3.046067	2.39E-09	6.9E-06
ADM	3.963986533	1.00703E-24	7.68616E-24	1.135633	2.02E-06	0.000777
FST	3.714380164	1.49909E-10	6.73409E-10	1.602592	6.56E-07	0.000345
BHLHB2	3.660332346	2.41494E-22	1.69131E-21	1.779734	1.77E-06	0.000722
CYR61	3.460940803	4.58412E-23	3.29964E-22	3.036142	5.96E-09	1.22E-05
ADRB2	3.044026009	2.97859E-16	1.67296E-15	2.437377	1.65E-07	0.000114
DUSP5	2.847290007	5.34772E-20	3.44196E-19	1.460276	1.57E-06	0.000673
HIST2H2AA4	2.724052988	1.01847E-13	5.18591E-13	1.200566	7.65E-07	0.00039
EFNB2	2.709520902	2.83635E-28	2.49632E-27	1.215793	3.36E-05	0.006438
HIST2H2AA3	2.675518726	4.78866E-12	2.2851E-11	1.228328	4.5E-06	0.001351
DNAJB1	2.57353698	4.89735E-39	7.0148E-38	1.344019	9.32E-06	0.002269
EGR1	2.560813139	2.52611E-12	1.21651E-11	3.634255	7.14E-10	3.83E-06
GADD45A	2.44572528	4.25018E-19	2.64854E-18	1.09097	5.53E-06	0.001573
CXCL1	2.252951442	1.28374E-05	4.36089E-05	2.16827	1.07E-07	9.02E-05
HAS3	2.187433524	4.91931E-08	1.94021E-07	1.039821	5.13E-05	0.008674
HERPUD1	2.134388824	1.48123E-20	9.71393E-20	1.002377	9.02E-06	0.002235
CCNL1	2.060797924	6.12652E-28	5.31988E-27	1.531865	2.08E-06	0.000777
MIR1974	2.053170518	1.13558E-10	5.13531E-10	2.579296	0.000258	0.025853
LOC100008589	1.967748936	7.25901E-12	3.44471E-11	1.843204	6.63E-07	0.000345
CCL20	1.794241775	0.000524911	0.001547247	1.158155	0.000195	0.021842
LOC100132564	1.732251157	3.18461E-07	1.19963E-06	1.95103	1.38E-07	0.000104
CD14	1.700324014	1.14268E-09	4.89735E-09	1.11706	7.12E-06	0.001878
ENC1	1.653764421	7.67012E-18	4.58172E-17	1.129762	0.000124	0.016402
CD59	1.586550226	3.87304E-23	2.79251E-22	1.081346	9.92E-05	0.01401
ATF3	1.475421346	4.40731E-09	1.83767E-08	4.302127	5.42E-10	3.76E-06
IFRD1	1.386579824	2.63634E-15	1.42823E-14	1.122456	1.43E-06	0.000627
IRS2	1.354135019	1.92067E-09	8.1541E-09	1.922806	4.6E-06	0.001363
BCL6	1.340384462	6.01622E-09	2.49152E-08	1.873812	4.7E-07	0.000263
GADD45B	1.263438308	5.43368E-14	2.79334E-13	2.55823	3E-08	2.98E-05
C13ORF15	1.258333669	2.41386E-05	8.03098E-05	2.651124	9.9E-08	8.61E-05
HIST3H2A	1.198692676	3.1811E-11	1.47308E-10	1.090508	5.43E-05	0.008964
IL6	1.193940996	0.000902549	0.00259803	1.953634	1.61E-06	0.000681
FOS	1.148625326	0.002655586	0.007273109	7.523147	6.92E-11	2.25E-06
CXCL2	1.13684321	2.58403E-05	8.5747E-05	3.846698	5.33E-08	5E-05
DDIT3	1.060898997	9.09349E-08	3.52906E-07	3.373708	3.22E-09	8.59E-06
CLK1	1.016256653	6.66236E-13	3.27574E-12	1.235296	3.43E-05	0.006521

Subsequently, we conducted a functional enrichment analysis for these 38 genes. GO analysis included 28 cSCC-PDT-related genes, which were mainly enriched in biological processes (BP), such as negative regulation of kinase activity, negative regulation of phosphorylation, and negative regulation of protein kinase activity ( ). Interestingly, DUSP6 participated in five different GO pathways ( ). KEGG analysis included 22 cSCC-PDT-related genes, which were significantly enriched in the MAPK signaling pathway ( ). To further evaluate the regulation of cSCC-PDT-related genes on the MAPK pathway, the results relevant to MAPK of GO analysis are highlighted in , indicating that cSCC-PDT-related genes are engaged in the MAPK pathway.

Figure 2

GO and KEGG enrichment analysis of 38 cSCC-PDT-related genes. (A) Radar plot of GO and KEGG analysis. (B) A chord gram showing the subordination between genes and the GO terms. (C) The list elaborates on the results of GO and KEGG analysis. (D) GO terms relevant to the MAPK pathway are highlighted.

Establishment and validation of diagnosis model

The Lasso regression was applied to establish a cSCC diagnostic model to further select key cSCC-PDT-related genes in GSE98767. Eventually, four of the 38 cSCC-PDT-related genes were included in the diagnostic model (). The ROC curve was used to assess the diagnostic utility of the diagnostic score and the expression level of DUSP6, EFNB2, DNAJB1, and CCNL1 in GSE98767 (cSCC and normal) and GSE45216 (cSCC and AK), and the boundary value with maximal ACU was selected as the cutoff (). In the training set, the diagnostic model had remarkable diagnostic utility (AUC = 1.000), while in the validation set, the overall diagnostic utility of the model was also noteworthy (AUC = 0.760).

DUSP6, EFNB2, DNAJB1, and CCNL1 were correlated with MAPK and cell cycle

The correlations between DUSP6, EFNB2, DNAJB1, CCNL1, and GOBP MAPK cascade gene sets were calculated by Spearman analysis. display a correlative coefficient between four cSCC-PDT genes and significant MAPK genes of the GOBP_MAPK_cascade in GSE98767 () and GSE45216 (). EFNB2 and DNAJB1 were reversely correlated with MAP3K6/MEKK6, MAPK3/ERK1, MAPK13/p38δ, and MAPK14/p38α. The PPI network exhibited the protein interactions of the above genes (). We then estimated the correlation between the four cSCC-PDT genes and the MAPK pathway promoter. The expression level of four cSCC-PDT genes was negatively correlated with some general MAPK pathway promoters, MRK, ERK, JNK, and p38 (11) in cSCC samples ().

Figure 3

Correlation analysis between DUSP6, EFNB2, DNAJB1, CCNL1, and MAPK pathway. (A) A heatmap of correlative coefficients between four cSCC-PDT genes and significant MAPK genes in GSE98767. (B) A heatmap of correlative coefficients between four cSCC-PDT genes and significant MAPK genes in GSE45216. (C) The PPI network of four cSCC-PDT genes and significant MAPK genes. (D–K) Spearman analysis of four cSCC-PDT genes and MAPK pathway promoter.

Previous studies have indicated that MAPK activation stimulates cell cycle transitions (14, 15). Here, we conducted ssGSEA on the cSCC sample of GES98767 to estimate the relationship between MAPK activation level and cell cycle. The normal group had a higher MAPK activation score, differentiation-enriched score, and lower cell cycle enriched score, angiogenesis-enriched score than the cSCC group (). Different MAPK promoters induce completely different functions, as JNK and p38 mediate cell apoptosis and downregulate in many tumors (27) whereas hyperactivation of ERK plays a crucial role in tumor proliferation (28). We wonder if the greater decline in JNK and p38 contributes to a lower MAPK activation score in cSCC. Consequently, the correlation between four cSCC-PDT-related genes and cell cycle genes was computed by spearman analysis in which DUSP6, EFNB2, DNAJB1, and CCNL1 were significantly correlated with significant cell cycle genes of GOBP_cell_cycle_process (). To explore how these genes regulate the MAPK signal pathway, we searched relevant signal maps from the KEGG database (https://www.genome.jp/kegg/), where DUSP6 inhibited ERK through dephosphorylation, thus inhibiting the EGF-EGFR-Ras-ERK signal pathway and preventing cell proliferation and differentiation ( ).

Figure 4

The four cSCC-PDT-related genes were correlated with the cell cycle gene. (A) The result of ssGSEA analysis of GSE98767 sorted by MAPK activation level. (B) Correlation coefficient matrix of four cSCC-PDT-related genes and cell cycle genes.

The four SCC-PDT related genes were potential PDT targets regulating the cell cycle by adjusting MAPK

To verify the relationship among the four SCC-PDT related genes, MAPK pathways, and cell cycle, we distinguished RAS, JNK, ERK, and P38 pathways in ssGESA analysis. The pathway annotation is shown in . At the pathway level of GSE98767, SCC samples had higher activity in the cell cycle, JNK_3 (negative regulation of JNK) than normal cells, while MAPK Activation Score, ERK and p38 were lower ( ). In the PDT-processed dataset, GSE84758, LC50 PDT-protocol treatment significantly increased the activity of MAPK Activation Score, RAS, JNK, ERK, and p38 (JNK, ERK, and p38 pathways included positive and negative regulatory pathways), as the LC90 PDT-protocol treatment reduced it ( ). The sample size of GSE84758 was too small to give a definitive inference, but enough to indicate some interesting phenomena.

Figure 5

Validation of the relationship among the four SCC-PDT related genes, the MAPK pathway, and the cell cycle. (A) A heatmap of ssGESA analysis on GSE98767, which is ordered by the expression level of DUSP6. (B) A heatmap of ssGESA analysis on GSE84758, which is ordered by the expression level of DUSP6. The treatment of A431 cells was annotated at the top. (C) The tSNA plot displayed the expression level of CCNL1, DNAJB1, DUSP6, and EFNB2. (D) PI staining and flow cytometry were performed to estimate the A431 cell cycle phase distribution. Recombinant EGF sc-4552 (100 ng/ml), ERK activator Ceramide C6 (10 µM), and JNK activator anisomycin (500 µM) were used to treat A431 cells before PDT treatment (2.8 J/cm2, 635 nm). The bar chart shows the composition ratio of cells in the G1 phase, S phase, and G2/M phase.

After that, we validated the above relationship at the single-cell level. The expression level of four cSCC-PDT-related genes in 5,902 HNSCC single cells of GSE103322 was displayed in the tSNE plot, and CCNL1 and DNAJB1 were highly expressed in most of the cells ( ). The gene expression level and cell cycle function state were comprehensively ranked, and the top five single cells with high cell cycle function states are shown in the radar plot ( ). To elucidate the effect of the MAPK pathway on cell cycles after PDT, A431 cells were incubated with 100 ng/ml recombinant EGF, 10 µM Ceramide C6 as the ERK overexpression group, and 500 nM anisomycin as the JNK overexpression group, respectively, before PDT treatment. The number of cells in each cell cycle was counted by flow cytometry. The results showed that more cells stagnated in the G1 phase in the JNK overexpression group compared to the control group ( ). Incidentally, ERK overexpression was supposed to promote cell proliferation and differentiation, producing effects similar to those of EGF treatment (29, 30). In this study, A431 cells with ERK overexpression performed similarly to cells with JNK overexpression after PDT treatment, resulting in inhibition of the cell cycle. In fact, different types of photosensitizers and cell lines activated MAPK pathways various. As EtNBSe is a newly developed photosensitizer, more pieces of evidence are needed to explore this phenomenon and its mechanism.

Discussion

The ascending incidence rates of cSCC increase the healthcare burden (31). PDT is recommended for treating cSCC because of its high tissue selectivity, good cosmetic outcome, comparatively less invasiveness, low recurrence rate, and, especially, safety for patients who cannot tolerate surgery. As for treating invasive cSCC, although many studies have demonstrated the efficacy of PDT (3, 32, 33), it has not been regarded as the first-line solution in the European guidelines (2). PDT still requires some modifications to be effectively used for invasive cSCC due to its inability to reach an adequate deep layer of skin. The combination of PDT and adjuvant therapy might expand the clinical application of PDT. Mechanism-based studies can provide evidence for the security and effectiveness of this combination. Previous studies illustrated that the MAPK pathway regulated cell proliferation, invasion, and differentiation in cSCC. PDT was proved to suppress the ERK pathway and activate the JNK pathway, and inhibit cell cycles. Nonetheless, there are rare reports about the comprehensive mechanism of PDT for cSCC treatment. Here we explored whether PDT inhibition of cSCC proliferation was related to MAPK signaling and its regulation of the G1 to G2/M phase.

In our study, we found 38 cSCC-PDT-related genes through differential analysis and intersection of differential gene sets. Then, we conducted GO/KEGG functional enrichment analysis for these 38 genes. The results indicated that most cSCC-PDT-related genes are enriched in BP, such as negative regulation of kinase activity, phosphorylation, and protein kinase activity. Many cSCC-PDT-related genes engaged in the MAPK pathway, which was consistent with previous studies (20, 34). Activation and inhibition of the MAPK pathway in cSCC are mainly affected by two aspects. On the one hand, the Ras/ERK signaling which is correlated to cancer initiation and development is upregulated in cSCC. Pierceall et al. found activating mutations and amplification of the Ras oncogene, which was an upstream activator of the Raf/Mek/Erk1/Erk2 kinase pathway, a typical MAPK pathway, in cSCC (35). Subsequent research determined that activating mutations in Ras could promote cSCC formation (36). Ratushny et al. clarified that MEK/ERK was one of the pro-oncogenic signaling pathways (37). On the other hand, JNK/p38 is downregulated so that the apoptosis of tumor cells is inhibited. Previous studies also illustrated these results (12, 16). In our study, MAPK activation was inhibited in tumor cells, probably because the downregulation of JNK was greater than the upregulation of ERK. As for the effect of PDT on the MAPK pathway, recent studies have reported that ZnPc-PDT (zinc phthalocyanine-PDT) could downregulate the expression and activity of MAPK (38) and PDT downregulated MAPK1 and MAPK3 (20). The above results indicated that the upregulation of the MAPK pathway might induce the occurrence and progression of cSCC, and PDT could suppress MAPK signaling and exert therapeutic effects. Additionally, it suggested that cSCC-PDT-related genes are probably engaged in the MAPK pathway. Previous research indicated that M-PDT could restrain MAPK activation by activating JNK (16). Furthermore, our flow cytometry results revealed that the mechanism of PDT inhibiting the proliferation of cSCC was that the upregulated JNK blocked the cell cycle transition from G1 to G2/M.

Next, we applied Lasso regression and established a diagnostic model based on four key cSCC-PDT-related genes, including DUSP6, EFNB2, DNAJB1, and CCNL1. After that, we preliminarily verified that the expression levels of four cSCC-PDT-related genes were negatively correlative with some MAPK pathway promoters by Spearman analysis. By ssGSEA analysis, we found that high MAPK activation was related to cell cycle transition. DUSP6 is a kind of dual-specificity mitogen-activated protein kinase phosphatase (MKP or DUSP), which is an established negative regulator of the MAPK pathway in mammalian cells (39). Groom et al. determined that DUSP6 preferentially inhibits ERK (40) and that the reduction of DUSP6 caused by ROS was correlated with high ERK activity (41). Additionally, the heat map showed high expression of DUSP6 in tumor cells compared with normal cells in , because DUSP6 was upregulated in response to elevated ERK signaling in cSCC with Ras activating mutation, where it was speculated to restrain ERK signaling (42, 43). Kim et al. found that EFNB2 inhibited vascular endothelial growth factor (VEGF)- and angiopoietin-1 (Ang1)-induced Ras/MAPK activities (44). In their study, Muller et al. observed a high prevalence of CCNL1 gene overexpression in HNSCCs and indicated that CCNL1 could be a regulator of cell-cycle transition (45). The function of DNAJB1 has rarely been reported. Furthermore, the single-cell analysis in our study also proved that CCNL1 and DNAJB1 were highly expressed in HNSCC.

Taken together, we speculated that increasing four cSCC-PDT-related genes inhibited the MAPK pathway in cSCC. Concretely, DUSP6 probably inhibited ERK5 MAPK signaling and EFNB2 restrained ERK1/2 and p38 MAPK signaling. DNAJB1 and CCNL1 were potential genes that repress JNK MAPK signaling. We innovatively established a connection between PDT, MAPK, and the cell cycle through bioinformatics analysis, proposing that PDT could repress the G1 to G2/M phase and inhibit tumor growth by regulating the MAPK pathway in cSCC. The results of flow cytometry were consistent with the abovementioned findings. There are still inescapable limitations to our research. Though we have screened for cSCC-PDT-related genes with bioinformatics analysis, attempting to clarify the roles of the MAPK pathway for treating PDT for cSCC, more in vivo and in vitro experiments are needed for the direct biological evidence to elucidate the exact mechanism of action. Meanwhile, we found a potential therapeutic target of PDT, which may be related to MAPK activation and cell cycle alternation. However, the definite effect of these four genes on cell cycles requires more exploration to expound.

Conclusion

CCNL1, DNAJB1, DUSP6, and EFNB2 were upregulated in cSCC or LC50 PDT-protocol treated A431 cells, which could be used as diagnostic markers for cSCC. Based on our results, we speculated that increased CCNL1, DNAJB1, DUSP6, and EFNB2 inhibited the MAPK pathway in cSCC. Despite a lower MAPK activation score, SCC appeared to have higher activity in the cell cycle. Strikingly, LC50 PDT-protocol treated A431 cell showed high levels of CCNL1, DNAJB1, DUSP6, and EFNB2 and extensive MAPK activation. Generally, we innovatively found a connection between PDT, MAPK and the cell cycle via CCNL1, DNAJB1, DUSP6, and EFNB2. In vitro, we discovered that PDT treatment suppressed the G1 to G2/M phase in JNK overexpressed A431 cells. In other words, CCNL1, DNAJB1, DUSP6, and EFNB2 are potential PDT target genes in cSCC treatment.

Data availability statement

The original contributions presented in the study are included in the article/ . Further inquiries can be directed to the corresponding authors.

Author contributions

YZ and JD wrote the first draft of the manuscript. YL, HW, and DZ contributed to conception and design of the study. YZ, JD, and YL contributed equally to this work. All authors contributed to the article and approved the submitted version.

Funding

This work is supported by the Hunan Medical Technology Innovation Guidance Fund grant no. 2020SK53615, the Natural Science Foundation of Hunan Province grant no. 2020JJ8015, 2019JJ40468, and 2019JJ80031, and the Changsha Natural Science Foundation grant no. kq2001041.

Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.