Identification of key genes in cutaneous squamous cell carcinoma: a transcriptome sequencing and bioinformatics profiling study.

BACKGROUND: Long-term exposure to ultraviolet (UV) radiation can cause cutaneous squamous cell carcinoma (cSCC), which is one of the most common malignant cancers worldwide. Actinic keratosis (AK) is generally considered a precancerous lesion of cSCC. However, the pathogenesis and oncogenic processes of AK and cSCC remain elusive, especially in the context of photodamage.
METHODS: In this study, transcriptome sequencing was performed on AK, cSCC, normal sun-exposed skin (NES) tissues, and normal non-sun-exposed skin (NNS) from 24 individuals. Bioinformatics analysis to identify the differentially expressed genes (DEGs) of 4 groups, and potential key genes of cSCC were validated by real-time quantitative reverse transcription PCR (qRT-PCR).
RESULTS: A total of 46,930 genes were differentially expressed in the 4 groups, including 127 genes that were differentially expressed between NES and NNS, 420 DEGs in AK compared to NES, 1,658 DEGs in cSCC compared to NES, and 1,389 DEGs in cSCC compared to AK. The Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis suggested that the DEGs are involved in multiple pathways, including extracellular matrix (ECM)-receptor interaction, immune, inflammatory, microbial infection, and other related pathways. Finally, 5 new genes (HEPHL1, FBN2, SULF1, SULF2, and TCN1) were confirmed significantly upregulated in cSCC.
CONCLUSIONS: Using transcriptome sequencing and integrated bioinformatical analysis, we have identified key DEGs and pathways in cSCC, which could improve our understanding of the cause and underlying molecular events of AK and cSCC. HEPHL1, FBN2, SULF1, SULF2, and TCN1 may be novel potential biomarkers and therapeutic targets of cSCC. 2021 Annals of Translational Medicine. All rights reserved.

Introduction

Cutaneous squamous cell carcinoma (cSCC) is one of the most common malignant cancers worldwide and is characterized by aberrant proliferation of keratinocytes (1). cSCC represents 20% of all non-melanoma skin cancers and is a significant cause of mortality due to its ability to metastasize to any organ in the body and its increasing prevalence (2). Actinic keratosis (AK) is one of the most common skin diseases in the elderly population and is generally considered a precancerous lesion of cSCC (3,4). Many studies have found that AK and cSCC exhibit highly similar clinical, histological, and molecular characteristics (5), and it is estimated that 10% of AK will develop into cSCC in approximately 2 years (6). The most important risk factor for cSCC and AK is exposure to ultraviolet B (UVB) radiation (7). UV radiation is considered a “complete carcinogen” because it is not only a mutagen but also a nonspecific damaging agent and has properties of both a tumor initiator and a tumor promoter (8). UV can cause cellular DNA damage and genetic mutations due to its powerful genotoxic effects (9,10). In addition, UV can induce inflammatory responses (11,12), promote the proliferation of epidermal keratinocytes (13), and cause immune tolerance or immunosuppression (14).

Studies have detailed the multiple mechanisms involved in the incidence and development of cSCC, including critical gene mutations (e.g., TP53, NOTCH 1-2, CDKN2A, FOS, RAS, EGFR, MYC, and BRAF) (15-18), key signaling pathway changes (e.g., PI3K/Akt, TGF-β, MAPK, NF-κB, and WNT) (19-22), immunodeficiency (23), microenvironment changes (24), epithelial-mesenchymal transition (EMT) (25), altered expression of adhesion molecules, reorganization of cytoskeletal proteins, and production of extracellular matrix (ECM)-degrading enzymes (26). However, skin cancer is a complex process, and the precise mechanisms involved in cSCC remain to be fully elucidated. Given the systemic metastatic and recurrent nature of cSCC (27), and the increasing incidence of AK and cSCC, more research is needed to identify key molecules involved in the transition of AK to cSCC, and the invasion and metastasis of cSCC. Furthermore, previous studies (28,29) examining the molecular pathogenesis of cSCC has primarily focused on fair-skinned populations in Europe, with few studies in Asian populations. Therefore, in this current study, transcriptome sequencing was conducted on 24 individuals from the plateau region of Yunnan, China, where UV causes a high incidence of AK and cSCC. Transcriptome sequencing and bioinformatics analysis can efficiently evaluate whole processes in one tissue in a global manner. Four types of tissues were examined, including AK, cSCC, normal sun-exposed skin (NES), and normal non-sun-exposed skin (NNS), and the differentially expressed genes (DEGs) were identified and further validated by real-time quantitative reverse transcription PCR (qRT-PCR). The results from this study provide an improved understanding of the occurrence and development of cSCC, and may provide helpful insights into early clinical defense mechanisms and effective treatment strategies for the management of AK and cSCC.

We present the following article in accordance with the MDAR reporting checklist (available at https://dx.doi.org/10.21037/atm-21-3915).

Methods

Tissue samples

Tissue samples of NNS (n=6), NES (n=6), AK (n=6), and cSCC (n=6) were obtained from patients during surgical excision at the Dermatology Department of the First Affiliated Hospital of Kunming Medical University (Yunnan, China). All AK and cSCC samples were derived from the face, while NES and NNS samples were derived from the face and trunk of healthy participants who were undergoing nevus resection surgery. Half of each sample was immediately snap-frozen in liquid nitrogen and stored at −80 °C for RNA extraction and gene expression analysis. The other half was fixed in formalin, embedded in paraffin, and used for histopathological diagnosis by two independent pathologists. Patients were aged 72.22±4.66 years (range, 65–78 years old). There were no significant differences in gender or age distribution among the four groups, and none of these patients had received any form of treatment prior to surgery. This study protocol was approved by the Ethics Committee of the First Affiliated Hospital of Kunming Medical University [Approval Number (2020)-L-29], and written informed consent was obtained from all patients. All procedures performed in this study involving human participants were in accordance with the Declaration of Helsinki (as revised in 2013).

RNA sequencing (RNA-Seq)

Total RNA extraction was performed using TRIzol reagent (TaKaRa Bio Inc., Kusatsu, Japan) following the manufacturer’s protocol, and RNA degradation and contamination were assessed by 1% agarose gel electrophoresis. Sequencing libraries were generated using the NEBNext®UltraTM RNA Library Prep Kit for Illumina® (NEB, USA), and library fragments were purified using the AMPure XP system (Beckman Coulter, Beverly, USA) to select cDNA fragments of the right length for PCR. Products were purified (AMPure XP system) and library quality was assessed on an Agilent Bioanalyzer 2100 system and sequenced on the Illumina HiSeq 4000 platform, and 150 bp paired-end reads were generated.

Bioinformatics analysis

The high-quality clean paired-end reads were aligned to the reference genome using Hisat2 software (http://ccb.jhu.edu/software/hisat2). The total number of mapped reads was determined, and the FPKMs (fragments per kilobase per million mapped reads) were calculated for each annotated gene. The R software package DESeq2 was employed to capture the DEGs and calculate the fold changes for each gene. Genes with|log2Fold change|>2, false discovery rate (FDR) <0.05, and Padj <0.05 were defined as DEGs and utilized for further analysis. In addition, the EMT-related genes (ERGs) and autophagy-related genes (ARGs) were obtained from public databases (http://dbemt.bioinfo-minzhao.org/index.html; http://www.tanpaku.org/autophagy/list/GeneList.html) to identify differentially expressed ERGs (DEERGs) and ARGs (DEARGs).

Function and pathway enrichment analysis

To analyze the identified DEGs at the functional level, significant Gene Ontology (GO, http://www.geneontology.org) biological process terms and Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg) were performed using the Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov) with the thresholds of Padj <0.05 and FDR <0.05. To identify pathways of intersecting genes, KEGG enrichment analysis was implemented. ClusterProfiler was applied to visualize the results.

Modules from the protein-protein interaction (PPI) network

To evaluate the interactive relationships among DEGs, the DEGs were mapped to the Search Tool for the Retrieval of Interacting Genes (STRING) database (https://string-db.org). Interactive DEGs were selected to construct the PPI network (combined score ≥0.75) and visualized using Cytoscape. The Molecular Complex Detection (MCODE) plugin in Cytoscape was used to screen the modules of the PPI network with MCODE score >3 and number of nodes ≥2.

Real-time quantitative reverse transcription PCR (qRT-PCR)

qRT-PCR was performed on an independent set of biopsies (n=5 in each group). Five significant DEGs were selected for qRT-PCR analysis, including TCN1, HEPHL1, FBN2, SULF1, SULF2, COL10A1, AREG, and ACTA1. Primers were designed and are listed in . Total RNA was isolated using TRIzol reagent (Invitrogen; Thermo Fisher Scientific, Inc.) and reverse transcribed to cDNA using an mRNA Reverse Transcription Kit (Invitrogen; Thermo Fisher Scientific, Inc.). qRT-PCR was performed using a SYBR-Green kit (Invitrogen; Thermo Fisher Scientific, Inc.). The RNA expression level of target genes was evaluated using 2−∆∆Ct.

Table 1

The primers of the differentially expressed genes used for qRT-PCR

Gene	Forward primer 5'-3'	Reverse primer 5'-3'
FBN2	TGAGTGCAGCATCATTCCTGG	AGAAACACATGCCTGTTCTCTG
HEPHL1	GAAAAACTACACCTACGTCTGGC	GATGTGCGAATGGTACACCCA
SULF1	TGGCGAGAATGGCTTGGATTA	TAACGGGCCTATGGGGATACA
SULF2	GGCAGGTTTCAGAGGGACC	GAAGGCGTTGATGAAGTGCG
TCN1	ATTGTCAGATGTAAGCTCGGGA	GCCATTGTGTGCTTCCATATTTT
GAPDH	GGACCTGACCTGCCGTCTAG	GTAGCCCAGGATGCCCTTGA

Statistical analysis

Statistical analyses were performed using SPSS Statistics 24.0 (IBM SPSS Statistics for Windows; IBM Corp, NY, USA). The data are presented as the means ± standard deviations. Differences among groups were analyzed using one-way analysis of variance (ANOVA) for normally distributed data and the Kruskal-Wallis test for non-normally distributed data. P<0.05 was considered statistically significant.

Results

Identification of the differentially expressed genes

The DEGs in the cSCC, AK, NES, and NNS samples were screened (https://cdn.amegroups.cn/static/public/atm-21-3915-1.xlsx). According to the screening criteria, a total of 127 DEGs (63 upregulated and 64 downregulated) were identified when comparing NES and NNS (), 420 DEGs (166 upregulated and 254 downregulated) were identified when comparing AK and NES samples (), 1,658 DEGs (485 upregulated and 1,173 downregulated) were identified when comparing cSCC and NES (), and 1,389 DEGs (448 upregulated and 941 downregulated) were identified when assessing cSCC and AK samples (). To identify the differentially expressed genes that were co-expressed in cSCC and AK, the DEGs detected during the following comparisons were overlapped: AK vs. NES, cSCC vs. NES, and AK vs. cSCC. As a result, 19 co-expressed DEGs were found (), including KRT32, KRT35, KRT85, CXCL3, CXCL11, ZNF536, CHIT1, F5, ACTA1, MUC16, CADM2, CSF3, CLVS2, SFRP4, LRRTM1, IL6, GFRA3, LINC02577, and EXTL1.

Figure 1

Analysis of the differentially expressed genes in the NNS, NES, AK, and cSCC samples. (A) Heat maps (left) and volcanic maps (right) of DEGs in the NES vs. NNS comparison; (B) heat maps (left) and volcanic maps (right) of DEGs in the AK vs. NES comparison; (C) heat maps (left) and volcanic maps (right) of DEGs in the cSCC vs. NES comparison; (D) heat maps (left) and volcanic maps (right) of DEGs in the cSCC vs. AK comparison; (E) venn diagram of the DEGs in the different groups. NNS, normal non-sun-exposed skin; NES, normal sun-exposed skin; AK, actinic keratosis; cSCC, cutaneous squamous cell carcinoma; DEG, differentially expressed gene.

Functional and pathway enrichment analysis of DEGs

The top 10 significant terms from the GO analysis are illustrated in https://cdn.amegroups.cn/static/public/atm-21-3915-2.xlsx. In the biological process (BP) terms, GO analysis results revealed that the upregulated DEGs in the NES vs. NNS analysis were primarily enriched in cell chemotaxis, antimicrobial humoral response, granulocyte chemotaxis and migration, and humoral immune response; while the downregulated DEGs were primarily enriched in the pattern specification process, regionalization, anterior/posterior pattern specification, embryonic skeletal system development and morphogenesis of a branching epithelium. In the AK vs. NES analysis, the upregulated DEGs were enriched in the regulation of lymphocyte activation, humoral immune response, regulation of chronic inflammatory response, and defense response to virus; while downregulated DEGs were enriched in the regulation of the ERK1 and ERK2 cascade and cellular response to metal ion. In the cSCC vs. NES analysis, upregulated DEGs were primarily enriched in the complement activation classical pathway and humoral immune response; and the downregulated DEGs were primarily enriched in skin development, keratinocyte differentiation, and regulation of membrane potential. In the cSCC vs. AK analysis, upregulated DEGs were primarily enriched in cell chemotaxis, skin development, endodermal cell differentiation, and connective tissue development; while downregulated DEGs were primarily enriched in keratinization, fatty acid metabolic process, regulation of lymphocyte activation, and lymphocyte and mononuclear cell proliferation.

For GO cell component (CC), upregulated DEGs in the NES vs. NNS analysis were primarily enriched in the cytoplasmic vesicle lumen, cell membrane, and specific granule; while downregulated DEGs were primarily enriched in the Wnt signalosome. In the AK vs. NES analysis, upregulated DEGs were primarily enriched on the external side of the plasma membrane and immunoglobulin complex; and downregulated DEGs were primarily enriched in the transcription factor AP-1 complex and collagen-containing extracellular matrix. For the cSCC vs. NES analysis, upregulated DEGs were primarily enriched in the extracellular matrix component, complex of collagen trimers, and basement membrane; and downregulated DEGs were primarily enriched in the intermediate filament cytoskeleton and transmembrane transporter complex. In the cSCC vs. AK analysis, upregulated DEGs were primarily enriched in the collagen-containing extracellular matrix, extracellular matrix component, endoplasmic reticulum lumen, basement membrane, and collagen trimer; while downregulated DEGs were primarily enriched in the intermediate filament, intermediate filament cytoskeleton, and keratin filament.

KEGG pathway enrichment analysis was used to identify the canonical signaling pathways associated with the significant DEGs. In the NES vs. NNS comparison, 42 significant KEGG pathways were annotated (). The upregulated DEGs were primarily involved in the interleukin (IL)-17 signaling pathway, tumor necrosis factor (TNF) signaling pathway, and cytokine-cytokine receptor interaction; while downregulated DEGs were primarily involved in signaling pathways regulating the pluripotency of stem cells, mTOR signaling pathway, Hippo signaling pathway, and Wnt signaling pathway. In the AK vs. NES comparison, 69 significant KEGG pathways were annotated (). The upregulated DEGs were primarily involved in the chemokine signaling pathway, viral protein interaction with cytokine and cytokine receptor, cytokine-cytokine receptor interaction, and Toll-like receptor signaling pathway; while downregulated DEGs were primarily involved in the IL-17 signaling pathway, TNF signaling pathway, and NOD-like receptor signaling pathway. In the cSCC vs. NES comparison, 60 significant KEGG pathways were annotated (). The upregulated DEGs were primarily involved in ECM-receptor interaction, human papillomavirus infection, PI3K-Akt signaling pathway, focal adhesion, mucin type O-glycan biosynthesis, regulating pluripotency of stem cells, and Wnt signaling pathway; while the downregulated DEGs were primarily involved in cell adhesion molecules (CAMs), PPAR signaling pathway, tyrosine metabolism, and cAMP signaling pathway. In the cSCC vs. AK analysis, 78 significant KEGG pathways were annotated (). The upregulated DEGs were mainly involved in ECM-receptor interaction, IL-17 signaling pathway, cytokine-cytokine receptor interaction, PI3K-Akt signaling pathway, focal adhesion, human papillomavirus infection, and TNF signaling pathway; and the downregulated DEGs were mainly associated with immunodeficiency and calcium signaling pathways.

Figure 2

Bubble plots of the top 10 enriched KEGG pathways of up- and downregulated differentially expressed genes. (A) Enriched pathways of the upregulated (left) and downregulated DEGs (right) identified in the NES and NNS samples; (B) enriched pathways of the upregulated (left) and downregulated DEGs (right) identified in the AK and NES samples; (C) enriched pathways of the upregulated (left) and downregulated DEGs (right) identified in the cSCC and NES samples; (D) enriched pathways of the upregulated (left) and downregulated DEGs (right) in the cSCC and AK samples. KEGG, Kyoto Encyclopedia of Genes and Genomes; NNS, normal non-sun-exposed skin; NES, normal sun-exposed skin; AK, actinic keratosis; cSCC, cutaneous squamous cell carcinoma; DEG, differentially expressed gene.

Protein-protein (PPI) network construction

Protein interactions among the DEGs were predicted using the STRING online database and Cytoscape software. A total of 37 nodes and 47 edges were involved in the PPI network of the DEGs identified in the NES vs. NNS analysis, 134 nodes and 333 edges were involved in the DEGs identified in the AK vs. NES analysis, 600 nodes and 2,618 edges were involved in the DEGs identified in the cSCC vs. NES comparison, and 462 nodes and 1647 edges were involved in the DEGs identified in the cSCC vs. AK comparison (https://cdn.amegroups.cn/static/public/atm-21-3915-3.xlsx). The results revealed that genes in the keratin family and the chemotactic cytokine family, including KRT6A, KRT6B, KRT6C, KRT16, KRT17, KRT28, KRT31, KRT32, KRT35, KRT79, KRT83, KRT85 and CXCL3, CXCL5, CXCL8, and CXCL11, were the most significantly altered genes in NES, AK, and cSCC tissues compared to NNS tissues.

Accumulating evidence has shown that EMT and autophagy play critical roles in cSCC (30,31). To understand the EMT and autophagy-related genes involved in the occurrence of cSCC, a total of 1,011 ERGs and 742 ARGs were screened from the RNA-Seq data and public databases and PPI maps were constructed for the significant DEERGs and DEARGs in cSCC compared to NES. The final PPI network complex contained 59 nodes and 141 edges ().

Figure 3

Protein-protein interaction network of DEERGs and DEARGs. All points represent differentially expressed mRNAs. The ellipse represents DEERGs, the rectangle represents DEARGs, and the parallelogram represents both DEARGs and DEERGs. Red represents upregulated DEGs, and green represents downregulated DEGs. DEERG, differentially expressed epithelial-mesenchymal transition-related gene; DEARG, differentially expressed autophagy-related gene.

Validation of significant differentially expressed genes using qRT-PCR

To validate the accuracy of the RNA sequencing results, qRT-PCR was performed on 5 significant DEGs in cSCC tissues, AK tissues, and NES tissues from 5 different patients. As shown in , the expression of HEPHL1, FBN2, SULF1, SULF2, and TCN1 was significantly upregulated in cSCC samples compared to NES and AK tissues, which was consistent with the RNA-Seq results.

Figure 4

Expression of HEPHL1, FBN2, TCN1, SULF1, and SULF2 as assessed by qRT-PCR in normal skin samples, AK samples, and cSCC samples. ns, P≥0.05; *P<0.05; **P<0.01; ***P<0.001. AK, actinic keratosis; cSCC, cutaneous squamous cell carcinoma.

Discussion

Although previous profiling experiments in cSCC have identified a large number of DEGs, the field has been hampered by a lack of consensus among these studies (3,32,33). Therefore, more in-depth studies are needed. Even as new sequencing technologies continue to evolve, transcriptome sequencing techniques still play a major role in identifying cancer-specific genes. In this study, next-generation sequencing technology was used to analyze the DEGs among NNS, NES, AK, and cSCC tissues from patients in Yunnan, China. A group of 5 DEGs that were highly associated with AK and cSCC tissues was identified and validated.

The gene expression patterns or the transcriptomes of cancer patients are greatly impacted by multiple factors, including infected microbiomes (34,35) and the severity of the disease (36), and thus, there can be considerable interpatient variation. Indeed, this study demonstrated that the expression levels of DEGs were significantly different between cSCC patients. Furthermore, the statistical results of all the DEGs in our study showed that each patient has a unique pattern of gene expression. The expression patterns of many DEGs in the transcriptome were not stable, which may be due to the heterogeneity of the tumor, the development of the disease, or the genetic background of the individual (37,38). This is a hindrance for researchers attempting to identify cancer driver genes for cSCC. Therefore, individual differences must be considered when conducting subsequent studies. Single-cell transcriptome sequencing (scRNA-seq) is a new technology for high-throughput sequencing analysis at the level of the single cell and may be a promising tool for exploring the heterogeneity and development of skin cancer (39).

By analyzing the overall changes in gene expression, the degree of DEG expression was most evident in AK and cSCC tissues, whereas normal skin appeared to have minimal alterations in gene expression. Based on human samples obtained from non-sun-exposed and sun-exposed skin, a small number of significant DEGs was identified. Among the top 20 significant DEGs, the S100 gene family (S100A8, S100A9, S100A12) and the SAA gene family (SAA1, SAA2) were significantly upregulated, while the HOX gene family (HOXA10, HOXB5, HOXB7, HOXA-AS3) was significantly downregulated. The S100 gene family plays an important role in the response to environmental triggers and cellular damage by regulating cellular differentiation, proliferation, apoptosis, inflammation, calcium homeostasis, cytoskeleton, and microbial resistance (40). S100 family proteins are mainly expressed in epidermal and skin tumors (41), among which S100A8 and S100A9 are considered to be key pro-inflammatory factors (42). S100A8 and S100A9 are secreted by psoriatic keratinocytes, and then send signals to adjacent keratinocytes through positive feedback to produce pro-inflammatory and pro-angiogenic cytokines, which exacerbate psoriasis skin lesions without affecting proliferation of keratinocytes. S100A8, S100A9, and S100A12 are positively correlated with cSCC tumorigenesis and highly expressed in cSCC. They can play an important role in the invasiveness of cSCC by enhancing in vitro proliferation and migration (41,43), which is consistent with our findings. The protein expression of the SAA gene is induced during acute and chronic inflammation or in response to tissue injury (44). Currently, SAA has been used as a predictor of cancer risk in various human malignancies. SAA is produced by the liver and enters the systemic circulation in response to the stimulation of inflammatory cytokines (such as IL-6) and participates in cholesterol transport, extracellular matrix degradation and inflammatory cell recruitment. The increase in SAA may be caused by too many inflammatory cytokines, especially IL-6. At present, not only has SAA highly expressed in a variety of cancer tissues, but also SAA gene families in AK and cSCC tissues are also significantly up-regulated in this experiment, indicating that elevated SAA levels may also be the main product of cancer (45,46). HOX genes encode a highly evolutionarily conserved transcription factor family known to be key regulators of embryonic development and various cellular processes, including cell growth, differentiation, apoptosis, motility, and angiogenesis (47). It can act as a tumor promoter or suppressor. Other DEGs, such as MMP3, CXCL5, CXCL8, CCL20, and CXCL1, were primarily enriched in inflammatory and immune-related pathways in AK and cSCC (48,49). Normal inflammation is self-limiting because the production of anti-inflammatory cytokines closely follows pro-inflammatory cytokines, while the subversion of cell death and/or repair processes occurs in chronically inflamed tissues, leading to DNA replication and loss of normal cell proliferation and growth control. Previous data demonstrated that the inflammatory response is a critical player in the tumor process, promoting the proliferation, survival, and migration of tumor cells (50). Therefore, the difference in the expression of DEG in AK and cSCC tissues suggests that inflammation may mediate the occurrence of the disease and occupy an important position.

The finding of similar DEGs between AK and cSCC tissues confirmed that AK is a precursor lesion of cSCC and indicated that they are genetically closely related. Nearly half of the DEGs identified in the AK vs. NES analysis overlapped with the DEGs identified in the cSCC vs. NES analysis. ECM-receptor interaction pathways were the most upregulated gene-enriched signaling pathways identified in the cSCC vs. NES analysis and the cSCC vs. AK analysis. The extracellular matrix (ECM) is reported to constitute the scaffold of tumor microenvironment and regulate cancer behavior (51). ECM-receptor interaction pathways were the most upregulated gene-enriched signaling pathways that play an important role in the process of tumor shedding, adhesion, degradation, movement and hyperplasia. The role of ECM in other cancers has been proved. ECM is upregulated in breast cancer tissue (52) and participates in the process of tumor invasion and metastasis in gastric cancer (53). Dysregulation of ECM remodeling caused by cancer cells promotes irreversible protein hydrolysis and cross-linking, which in turn affects cell signaling, microenvironmental cues, angiogenesis, and tissue biomechanics (54). Matrix metalloproteinases (MMPs) are an important example. MMPs are zinc-containing endopeptidases with an extensive range of substrate specificities. These enzymes are able to degrade various components of ECM proteins and are involved in angiogenesis, which promotes cancer cell growth and migration (55). In our results, increased expression of MMP13, MMP10, and MMP1 was observed in cSCC samples, which was consistent with previous reports (56,57). This indicated that MMP might affect the ECM-receptor interaction pathway by degrading various ECM proteins to promote the occurrence and development of tumors. Microorganisms exist on the surface of skin and participate in various physiological processes. All microorganisms living in the human body are an indispensable part of maintaining health and immune system (58). Recent studies demonstrate correlations between a particular composition of microbiome and a broad spectrum of diseases, including lung cancer, obesity, and even breast cancer (59,60). Microbiome-related pathways were also enriched in both AK and cSCC, especially the human papilloma virus (HPV), which has been shown to be positively correlated with keratinocyte cell proliferation, inflammation, and the immune response (34,61), suggesting that HPV is closely related to the occurrence of AK and cSCC and the transformation of AK to cSCC (62-64). In addition, a group of EMT- and autophagy-related genes involved in cSCC was identified, and most of the DEGs related to autophagy were downregulated. These results suggested that inhibition of autophagy and activation of EMT plays an important role in the occurrence and development of cSCC.

Moreover, 5 highly expressed DEGs were confirmed by qRT-PCR in the cSCC vs. NES analysis and the cSCC vs. AK analysis. Both hephaestin-like 1 (HEPHL1), which is thought to function as a ferroxidase and facilitates the transfer of iron to transferrin, and HMOX1, were highly expressed in cSCC. Only a few studies have focused on the role of HEPHL1 in tumors. HEPHL1 mutation has been detected in POEMS syndrome (65) and may also be related to the malignant transformation of hepatocyte L02 cells induced by microcystin-LR (66). Therefore, what role HEPHL1 plays in cSCC needs to be further studied. Similar to HEPHL1, FBN2, TCN1, SULF1, and SULF2 were also highly expressed in cSCC. Fibrillin 2 (FBN2) encodes a component of connective tissue microfibrils and forms the structural core of microfibrils. Previous studies have shown that methylation or silencing of FBN2 in tumor cells may play an important role in carcinogenesis, invasion, and metastasis (67). However, recent research suggested the opposite, wherein FBN2 stimulates tumor angiogenesis by influencing TGF-β sequestration by microfibrils, leading to a locally higher active TGF-β concentration in the tumor microenvironment which is involved in tumor progression, metastasis, EMT, and tumor angiogenesis (68). In lung cancer, cholangiocarcinoma, and metastatic colorectal cancer, FBN2 has been found to be associated with EMT, cell migration, and metastasis (69-71). In combination with our KEGG pathway enrichment analysis, we found that the ECM-receptor interaction pathway was significantly up-regulated and the expression of MMP protein increased. We speculate that FBN2 may be their inducing factor, and affect the tumor microenvironment by acting on them. Transcobalamin I (TCN1) encodes a member of the vitamin B12-binding protein family, is expressed in various tissues and secretions, and regulates cobalamin homeostasis. Studies have shown that TCN1 is associated with poor outcomes in malignancies, such as pancreatic and colon cancer (72,73). Sulfatase 1 (SULF1) and sulfatase 2 (SULF2) encode extracellular heparan sulfate endosulfatase, which can edit the sulfation status of heparan sulfate proteoglycans (HSPGs) on the outside of cells and regulate a number of critical signaling pathways, such as Wnt/β-catenin signaling (74). Recent studies have shown that SULF1 and SULF2 are involved in the occurrence and development of many tumors, including breast cancer, lung cancer, and prostate cancer, indicating their potential for cancer diagnosis and therapy (75-77). SULF1 is highly expressed in colon cancer, and an in vitro study confirmed that M2-like macrophages in colon cancer increased the amount of SULF1 produced by fibroblasts and were associated with poor prognosis (78). Graham et al. (79) found that the activation of SULF1/SULF2 not only promoted fetal growth regulation and injury-induced liver regeneration but also dysregulated tumor growth. Interestingly, SULF1 and SULF2 were significantly elevated in cSCC compared to AK and normal skin in our study, indicating that the SULF family plays an important role in cSCC and may mediate tumor growth.

Interestingly, the expression of ACTA1 in NES, AK, and cSCC tissues exhibited a downward trend as the disease progressed. ACTA1 is downregulated in head and neck squamous cell carcinoma and may serve as a potential diagnostic and prognostic marker of tumors (80). A similar result was found in colorectal cancer and pancreatic adenocarcinoma (81). A large amount of research shows that keratin is associated with a variety of diseases (82,83). In the current study, the keratin gene family, including KRT 85, KRT 32, and KRT 35, displayed the same downward trend as ACTA1. Actin and keratin are known to be components of the cytoskeleton (84) (a schematic diagram of the cytoskeleton is shown in ). The various components of the cytoskeleton are well orchestrated in normal cells, and mutations and abnormal expression of cytoskeletal and cytoskeletal-associated proteins play an important role in cancer cell migration, invasion, and metastasis. The actin family of proteins are highly conserved proteins that play a role in cell motility, structure, and integrity. Studies have shown that actin cytoskeletal rearrangement is associated with invasion and metastatic phenotypes of various cancer cells (85-87). Furthermore, upregulation of the keratin family KRT80 can promote cytoskeletal rearrangements at the leading edge and increase focal adhesion and cellular stiffening, collectively promoting cancer cell invasion that is associated with poor survival (88). High expression of KRT80 in gastric cancer and colorectal cancer is also involved in the malignant proliferation of tumors (89,90). A study of basal cell carcinoma (BCC) showed a close correlation between differentiation‐specific keratins and BCC (91). Our study found that keratinocyte differentiation and hyperproliferation-related genes, and stress markers (e.g., KRT6A, KRT6B, KRT6C, KRT14, and KRT16) were upregulated in cSCC. Some members of the keratin family (KRT32, KRT35, KRT85) were significantly and consistently downregulated in AK and cSCC compared to normal skin. This could suggest that certain types of keratin may play different roles in the progression of cancer, especially in the migration, invasion and metastasis of cancer cells.

Figure 5

Schematic diagram of the cytoskeleton.

CXCL11 is also highly expressed in cSCC, and studies have shown that it can bind to CXCR3 and CXCR7, which is associated with angiogenesis, invasiveness, and reduction of apoptosis in tumor cells (92,93). GFRA3 is a glycosylphosphatidylinositol (GPI)-linked cell surface receptor. Eftang et al. (94) found that high methylation levels of GFRA3 conferred a very unfavorable prognosis in gastric cancer. However, expression of GFRA3 was associated with the growth of cancer cells, tumor lymph node metastases, and higher clinical stages of breast cancer and non-small cell lung cancer (95,96). These studies suggested that GFRA3 may play different roles in various cancer types. In our study, GFRA3 was downregulated and may be related to antitumor properties, but the exact mechanisms remain to be further studied.

Conclusions

The molecular characterization of AK and cSCC is an important step towards identifying and understanding the dysregulated pathways and key drivers of cancer development and progression. This study identified some important genes that may play crucial roles in the occurrence and development of cSCC, and perhaps AK. To the best of our knowledge, this is the first study to examine the function of these genes in cSCC. Additional studies should be performed to further explore the pathogenic mechanisms of these genes.

The article’s supplementary files as