Get 20M+ Full-Text Papers For Less Than $1.50/day. Subscribe now for You or Your Team.

Learn More →

Identifying an lncRNA-Related ceRNA Network to Reveal Novel Targets for a Cutaneous Squamous Cell Carcinoma

Identifying an lncRNA-Related ceRNA Network to Reveal Novel Targets for a Cutaneous Squamous Cell... biology Article Identifying an lncRNA-Related ceRNA Network to Reveal Novel Targets for a Cutaneous Squamous Cell Carcinoma 1 1 1 , 2 , 3 3 1 Yaqin Xu , Yingying Dong , Yunhua Deng *, Qianrong Qi * , Mi Wu , Hongmei Liang , Qiuyun She and Qing Guo Department of Dermatology, Tongji Medical College, Huazhong University of Science and Technology, Wuhan 430030, China; xuyaqin@hust.edu.cn (Y.X.); yingyingdong@hust.edu.cn (Y.D.); qiuyunshe@hust.edu.cn (Q.S.) Department of Obstetrics & Gynecology, University of California, Irvine, CA 92697, USA Department of Immunology, School of Basic Medicine, Tongji Medical College, Huazhong University of Science and Technology, Wuhan 430030, China; wumi@hust.edu.cn (M.W.); lianghongmei@hust.edu.cn (H.L.); qingguo@hust.edu.cn (Q.G.) * Correspondence: yhdeng@tjh.tjmu.edu.cn (Y.D.); qianronq@hs.uci.edu (Q.Q.) Simple Summary: The exact functions and molecular mechanism of lncRNAs, acting as competitive endogenous RNAs in a cutaneous squamous cell carcinoma, remain unexplored. The present study was conducted to identify the differentially expressed lncRNAs and mRNAs and establish the lncRNA-related competing endogenous RNA networks associated with a cutaneous squamous cell carcinoma. A competing endogenous RNA network consisting of 137 miRNA-lncRNA and 221 miRNA-mRNA pairs was constructed. As for the functional analysis of the mRNAs in the network, a FoxO signaling pathway, an autophagy and cellular senescence were the top enrichment terms based on the Kyoto Encyclopedia of Genes and Genomes analysis. We identified five core Citation: Xu, Y.; Dong, Y.; Deng, Y.; mRNAs and built a core mRNA-associated competing endogenous RNA network. Finally, one Qi, Q.; Wu, M.; Liang, H.; She, Q.; lncRNA HLA-F-AS1 and three mRNAs named AGO4, E2F1 and CCND1 in the core mRNA-associated Guo, Q. Identifying an lncRNA- competing endogenous RNA network were validated with the same expression patterns. The core Related ceRNA Network to Reveal mRNAs and their associated lncRNAs may provide potential therapeutic targets for cutaneous Novel Targets for a Cutaneous squamous cell carcinomas. Squamous Cell Carcinoma. Biology 2021, 10, 432. https://doi.org/ Abstract: A cutaneous squamous cell carcinoma (cSCC) derived from keratinocytes is the second 10.3390/biology10050432 most common cause of non-melanoma skin cancer. The accumulation of the mutational burden of genes and cellular DNA damage caused by the risk factors (e.g., exposure to ultraviolet radiation) Academic Editor: Loredana Moro contribute to the aberrant proliferation of keratinocytes and the formation of a cSCC. A cSCC Received: 22 March 2021 encompasses a spectrum of diseases that range from recursor actinic keratosis (AK) and squamous Accepted: 9 May 2021 cell carcinoma (SCC) in situ (SCCIS) to invasive cSCCs and further metastatic SCCs. Emerging Published: 13 May 2021 evidence has revealed that lncRNAs are involved in the biological process of a cSCC. According to the ceRNA regulatory theory, lncRNAs act as natural miRNA sponges and interact with miRNA Publisher’s Note: MDPI stays neutral response elements, thereby regulating the mRNA expression of their down-stream targets. This with regard to jurisdictional claims in study was designed to search for the potential lncRNAs that may become potential therapeutic published maps and institutional affil- targets or biomarkers of a cSCC. Considering the spirit of the study to be adequately justified, we iations. collected microarray-based datasets of 19 cSCC tissues and 12 normal skin samples from the GEO database (GSE42677 and GSE45164). After screening the differentially expressed genes via a limma package, we identified 24 differentially expressed lncRNAs (DElncRNAs) and 3221 differentially expressed mRNAs (DEmRNAs). The miRcode, miRTarBase, miRDB and TargetScan databases were Copyright: © 2021 by the authors. used to predict miRNAs that could interact with DElncRNAs and DEmRNAs. A total of 137 miRNA- Licensee MDPI, Basel, Switzerland. lncRNA and 221 miRNA-mRNA pairs were retained in the ceRNA network, consisting of 31 miRNAs, This article is an open access article 11 DElncRNAs and 155 DEmRNAs. For the functional analysis, the top enriched biological process distributed under the terms and was enhancer sequence-specific DNA binding in Gene Ontology (GO) terms. The FoxO signaling conditions of the Creative Commons pathway, autophagy and cellular senescence were the top enrichment terms based on a Kyoto Attribution (CC BY) license (https:// Encyclopedia of Genes and Genomes (KEGG) analysis. The combination of a STRING tool and creativecommons.org/licenses/by/ Cytoscape software (plug-in MCODE) identified five core mRNAs and built a core mRNA-associated 4.0/). Biology 2021, 10, 432. https://doi.org/10.3390/biology10050432 https://www.mdpi.com/journal/biology Biology 2021, 10, 432 2 of 15 ceRNA network. The expression for five identified core mRNAs and their related nine lncRNAs was validated using the external dataset GSE7553. Finally, one lncRNA HLA-F-AS1 and three mRNAs named AGO4, E2F1 and CCND1 were validated with the same expression patterns. We speculate that lncRNA HLA-F-AS1 may sponge miR-17-5p or miR-20b-5p to regulate the expression of CCND1 and E2F1 in the cSCC. The present study may provide potential diagnostic and therapeutic targets for cSCC patients. Keywords: competing endogenous RNAs network; long non-coding RNAs; microRNAs; cutaneous squamous cell carcinoma 1. Introduction Cutaneous squamous cell carcinoma (cSCC), which originates from epidermal ker- atinocytes, ranks as the second most occurring skin cancer [1]. The occurrence of cSCC is due to many risk factors. Beside the most important risk factor, cumulative ultraviolet exposure, other risk factors include old age, immunosuppression, smoking and genetic factors [2]. Moreover, a cSCC commonly has high mutational burden (e.g., TP53, CDKN2A, Ras and NOTCH1). The accumulation of these mutations and cellular DNA damage caused by the risk factors contribute to the uncontrolled proliferation of keratinocytes and the formation of a cSCC and disease progression [3]. Consequently, a cSCC is a spectrum of diseases that commonly arises from precursor actinic keratosis (AK) and a squamous cell carcinoma (SCC) in situ (SCCIS) to an invasive cSCC to a metastatic SCC with a distant spread [4]. Based on the tumor stage, the optimal treatment approach for a cSCC is selected. The cSCC at the early stage is generally identified and can be cured via surgical resection and/or radiotherapy. For cases at the advanced stage (invasive and metastatic cSCCs), a single or a combined therapy, including chemotherapy, radiation therapy and epidermal growth factor receptor (EGFR) inhibitors, have been implemented with limited and unsat- isfactory efficacy with a long-term survival rate of less than 20% [4]. Thus, it is urgent to identify new biomarkers and targets that are more effective for the early diagnosis and treatment of a cSCC. Non-coding RNAs (ncRNAs) incapable of protein coding can be classified according to their length including long non-coding RNAs (lncRNAs) (>200 nucleotides) and small ncRNAs (<200 nucleotides) such as microRNAs (miRNAs), small nucleolar RNAs and transfer RNAs [5]. The dysregulation of the expression of miRNAs has also been elucidated to be implicated in the tumorigenesis of a cSCC. The upregulated miR-142-5p induced cancer stem cell-like properties in a cSCC via the direct regulation of PTEN [6]. MiR-27a directly bound to EGFR, leading to the inactivation of the NF-B pathway and the inhibition of cSCC cell growth [7]. LncRNAs are a wide and diverse class of single-stranded ncRNAs longer than 200 nu- cleotides lacking a protein-coding capacity [8]. LncRNA dysregulation has also been linked to the process of cell differentiation and growth and the pathogenesis of many types of cancers [9,10] including cSCCs [11,12]. LINC00520 targeted EFGR to inactivate the PI3K/Akt pathway and restrain the cSCC progress [11]. The lncRNA PICSAR dys- regulated the collagen and fibronectin expression, resulting in the dysregulation of the adhesion and migration of cSCC cells [12]. Emerging evidence demonstrates that activities of lncRNAs and miRNAs are related to each other in different complex mechanisms [13]. One of these mechanisms is the role of lncRNAs as competing endogenous RNAs (ceR- NAs). Specifically, lncRNAs can compete with mRNA to bind miRNA binding sites, thereby negatively regulating miRNAs and their target genes. In other words, lncRNAs act as a sponge for miRNAs abolishing their inhibitory action to target mRNAs [14]. The lncRNA-miRNA-mRNA regulatory network has been verified to be associated with the progression of several types of squamous cell carcinomas (SCCs) [15–17]. LncRNA nuclear paraspeckle assembly transcript 1 (NEAT1) sponged to miR-129, leading to the depression Biology 2021, 10, 432 3 of 15 of CTBP2 and regulated the esophageal squamous cell carcinoma (ESCC) cell viability and invasion [15]. LINC00511 induced the expression of LAMC2 via sponging to miR-765 and modulated a tongue squamous cell carcinoma progression (TSCC) [16]. Zhao et al. indicated that lncRNA HCP5 functioned as a molecular sponge to absorb miR-140-5p in an oral squamous cell carcinoma (OSCC) and upregulated the SOX4 gene expression [17]. As for cSCCs, lncRNA HOTAIR promoted the progression of a cSCC, exerting a function as a ceRNA through the miR-326/Prenylated Rab acceptor 1 domain family member 2 (PRAF2) axis [18]. In this paper, we compared total RNA expression data based on an array between cSCC and normal skin (NS) samples from the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database. Based on the differentially expressed genes (DEGs), we performed a series of down-stream analyses including the construction of an lncRNA-related ceRNA network, a functional enrichment analysis and a protein- protein interaction analysis. We theorize that lncRNA HLA-F-AS1 may act as the sponge of miR-17-5p and miR-20b-5p to regulate the expression of CCND1 and E2F1 in the cSCC. The expression of HLA-F-AS1 in the cSCC samples and its role mechanism in the biological process of a cSCC as a ceRNA required further research. They may be candidate targets with a promising diagnostic or therapeutic value with a cSCC. The remainder of this paper is organized as follows. Section 2 introduces the dataset collection and presents our designed data processing workflow and describes our down- stream analyses and validation methods. The experiment results are described in Section 3. The Discussion based on the experiment results is analyzed in Section 4. Conclusions are drawn in Section 5. 2. Materials and Methods 2.1. Dataset Collection from the GEO Database The expression profiles of lncRNAs and mRNAs between a cSCC and normal skin (NS) from healthy samples were downloaded from the GEO database (www.ncbi.nlm.nih.gov/ geo, accessed on 24 August 2020). The criterion for selecting the cSCC samples was that the clinical and histological diagnosis of biopsies obtained during surgery was a primary cSCC (SCCIS or invasive cSCC). The datasets were recruited with the standard based on the same microarray platform and containing more than three samples in each group. Finally, GSE42677 and GSE45164 dataset bases on the GPL571 platform were obtained including 20 cSCC and 13 NS samples. 2.2. Data Preprocessing and Screening for DElncRNAs and DEmRNAs Based on the annotation information from Affymetrix Human Genome U133A 2.0 Array platform, the probe sequences were quantified to the expression of lncRNAs and mRNAs after mapping to the human GRCh38 reference genome (http://useast.ensembl. org/info/data/ftp/index.html, accessed on 29 August 2020). We integrated all samples from the two datasets (GSE42677 and GSE45164) and applied ComBat [19] to remove the batch-effects of expression values (Figure S1a,b). A principal component analysis (PCA) [20] was applied to assess the independence of the available genes between 20 cSCC and 13 NS samples. According to the PCA result (Figure S1), we removed one cSCC and one control sample affecting the classification. Data of 19 cSCC patients and 12 NS were selected for a subsequent analysis (Figure 1a). Biology 2021, 10, x 4 of 16 Volcano plots were performed for the DEGs using the ggplot2 package [22] (version 3.3.2; Biology 2021, 10, 432 4 of 15 https://cran.r-project.org/web/packages/ggplot2/index.html, accessed on 29 August 2020) in the R software. Pheatmap package (version: 1.0.12; https://cran.r-project.org/web/pack- ages/pheatmap/, accessed 29 August 2020) in R was used to generate the hierarchical cluster heatmap representing the expression intensity and direction (Figure 1b). (a) Biology 2021, 10, x 5 of 16 (b) (c) Figure 1. Differentially expressed analysis between cSCC and NS tissues from GSE42677 and GSE45164 datasets. (a) The principal component analysis (PCA) of the selected 19 cSCC and 12 NS tissues; (b) volcano plots showing DEGs. The red Figure 1. Differentially expressed analysis between cSCC and NS tissues from GSE42677 and GSE45164 datasets. (a) The dots represent an upregulated expression and the blue dots represent downregulated DEGs; (c) cluster analysis (heatmap) between 19 cSCC and 12 NS tissues. The red stripes display upregulation while the blue stripes display the downregula- principal component analysis (PCA) of the selected 19 cSCC and 12 NS tissues; (b) volcano plots showing DEGs. The red tion of genes. DEGs: differential expression genes. dots represent an upregulated expression and the blue dots represent downregulated DEGs; (c) cluster analysis (heatmap) 2.3. Construction of the cSCC-Associated lncRNA-miRNA-mRNA Network between 19 cSCC and 12 NS tissues. The red stripes Due to displ the absence ay upr of egulation miRNA information while , the the miRco blue de data stripes base [23] display the downregulation (http://www.mircode.org/, accessed on 1 September 2020) was used to predict the targeted miRNAs of DElncRNAs. According to the miRTarBase [24] (http://mirtar- of genes. DEGs: differential expression genes. base.mbc.nctu.edu.tw/, accessed on 1 September 2020), miRDB [25] (http://www.mirdb.org/, accessed on 1 September 2020) and TargetScan [26] (http://www.targetscan.org/), the DEmRNA-miRNA interaction pairs were obtained. Only the overlapping of the interaction pairs supported by all databases was retained to create the miRNA-mRNA pairs. Subsequently, the correlation between DElncRNAs and DEmRNAs was analyzed by Pearson’s correlation coefficient (PCC) with the threshold of the PCC value r > 0.4 and p < 0.05 indicating a strong correlation. Finally, the remaining Biology 2021, 10, 432 5 of 15 A Limma package [21] (version 3.44.3; https://bioconductor.org/packages/release/ bioc/html/limma.html, accessed on 29 August 2020) in R was applied to identify the differentially expressed lncRNAs (DElncRNAs) and differentially expressed mRNAs (DEmRNAs) before a quantile normalization and base-2 logarithmic (log ) transforma- tion. The criteria were selected as a false discovery rate (FDR) < 0.05 and a fold change (FC) > 1.5. Volcano plots were performed for the DEGs using the ggplot2 package [22] (version 3.3.2; https://cran.r-project.org/web/packages/ggplot2/index.html, accessed on 29 August 2020) in the R software. Pheatmap package (version: 1.0.12; https://cran.r- project.org/web/packages/pheatmap/, accessed 29 August 2020) in R was used to gen- erate the hierarchical cluster heatmap representing the expression intensity and direction (Figure 1b). 2.3. Construction of the cSCC-Associated lncRNA-miRNA-mRNA Network Due to the absence of miRNA information, the miRcode database [23] (http://www. mircode.org/, accessed on 1 September 2020) was used to predict the targeted miR- NAs of DElncRNAs. According to the miRTarBase [24] (http://mirtarbase.mbc.nctu. edu.tw/, accessed on 1 September 2020), miRDB [25] (http://www.mirdb.org/, accessed on 1 September 2020) and TargetScan [26] (http://www.targetscan.org/), the DEmRNA- miRNA interaction pairs were obtained. Only the overlapping of the interaction pairs supported by all databases was retained to create the miRNA-mRNA pairs. Subsequently, the correlation between DElncRNAs and DEmRNAs was analyzed by Pearson’s correlation coefficient (PCC) with the threshold of the PCC value r > 0.4 and p < 0.05 indicating a strong correlation. Finally, the remaining miRNA-DElncRNA and miRNA-DEmRNA and their shared miRNAs were incorporated into the ceRNA network. Cytoscape [27] (version 3.6.1; http://cytoscape.org, accessed on 6 September 2020) was utilized to visualize the network. 2.4. Functional Enrichment Analysis The Gene Ontology (GO) functional enrichment analysis in the category Biologi- cal Processes (BP) and the KEGG pathway enrichment analysis of DEmRNAs in the ceRNA network were determined by the clusterProfiler package in R [28] (version 3.14.0; https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html, accessed on 12 September 2020) and R. The threshold was selected to be p < 0.05. 2.5. Establishment of the Protein-Protein Interaction (PPI) Network and the Core ceRNA Subnetwork The Search Tool for the Retrieval of Interacting Genes (STRING) [29] (http://string- db.org/, accessed on 12 September 2020) was utilized to analyze PPI information with the threshold of a combined score >0.4. The plug-in Molecular Complex Detection (MCODE; version:1.4.2; http://apps.cytoscape.org/apps/mcode, accessed on 12 September 2020) [30] in Cytoscape was utilized to identify functionally-related hub modules with a degree cut- off of more than 5, a node score cut-off of 0.2, a k-score of 2 and maximum depth of 100. The key genes and related lncRNAs and miRNAs were then mapped into the core ceRNA subnetwork. 2.6. Expression Validation in an External Dataset To validate the expression of the DElncRNAs and DEmRNAs in the subnetwork, the GSE7553 based on the GPL570 platform with 11 SCC and 4 NS samples was obtained from GEO. The inclusion criteria of the cSCC samples, the data processing of the identifying DElncRNAs and DEmRNAs of GSE7553 datasets and the cut-off criteria were consistent with those of the GSE42677 and GSE45164 datasets. If the relative expression level of genes in the core ceRNA network was significantly different and met the same expression patterns with the GSE42677 and GSE45164 datasets, we considered those genes validated. The comparison of the log transformed expression level of the DElncRNAs and DEmRNAs was performed by GraphPad Prism V8.3.0 software (GraphPad Software, Inc., San Diego, CA, USA) with an unpaired t-test and p < 0.05 indicating a significant difference. Biology 2021, 10, 432 6 of 15 3. Results 3.1. A Total Of 24 DElncRNAs and 3221 DEmRNAs Were Identified between cSCC and NS Based on the PCA results of all samples from the GSE42677 and GSE45164 datasets, one cSCC sample and one NS sample affecting the classification were excluded for the down- stream analyses (Figure S1c). The heatmap of the DEGs showed that the gene expression levels between the remaining 19 cSCC and 12 NS samples were completely distinguishable (Figure 1c). Compared with the NS samples, 24 DElncRNAs (18 over-expressed and 6 down-expressed) and 3221 DEmRNAs (1568 over-expressed and 1653 down-expressed) were identified with the threshold of FDR < 0.05 and |FC| > 1.5. 3.2. The cSCC-Associated ceRNA Network Establishment Before predicting the miRNAs interacting with DElncRNAs or DEmRNAs, the co- expression analysis by the PCC was applied to calculate the correlation between DElncRNA and DEmRNA. Based on the result of the PCC with a cut-off value r > 0.4 and p < 0.05, 155 DEmRNAs and 23 DElncRNAs were retained for a further analysis (Table S1). The re- maining DEmRNAs and DElncRNAs that shared the common miRNAs were incorporated Biology 2021, 10, x into the construction of the ceRNA network (Figure 2). Finally, a cSCC-associated 7 ceRNA of 16 network containing 137 miRNA-lncRNA (Table S2) and 221 miRNA-mRNA (Table S3) pairs was created composed of 31 miRNAs, 11 DElncRNAs (8 down-expressed and 3 over- expressed, Table S4) and 155 DEmRNAs (130 down-expressed and 25 over-expressed) and 3 over-expressed, Table S4) and 155 DEmRNAs (130 down-expressed and 25 over- (Tables S4 and S5). expressed) (Tables S4 and S5). Figure 2. The construction of the ceRNA network including 11 DElncRNAs, 31 predicted miRNAs and 155 DEmRNAs. The Figure 2. The construction of the ceRNA network including 11 DElncRNAs, 31 predicted miRNAs and 155 DEmRNAs. The circle in green, the rectangle in red and the triangle in blue represent DEmRNAs, miRNAs and DElncRNAs, respectively. circle in green, the rectangle in red and the triangle in blue represent DEmRNAs, miRNAs and DElncRNAs, respectively. DElncRNAs: differentially expressed lncRNAs; DEmRNAs: differentially expressed mRNAs; miRNAs: microRNAs. DElncRNAs: differentially expressed lncRNAs; DEmRNAs: differentially expressed mRNAs; miRNAs: microRNAs. 3.3. DEmRNAs of the ceRNA Network Were Associated with Enhancer-Binding, Protein Kinase 3.3. DEmRNAs of the ceRNA Network Were Associated with Enhancer-Binding, Protein Kinase Regulator Activity and Autophagy Regulator Activity and Autophagy To uncover the biological mechanisms associated with the development of a cSCC, a To uncover the biological mechanisms associated with the development of a cSCC, GO functional and a KEGG pathway enrichment analysis of the DEmRNAs in the ceRNA a GO functional and a KEGG pathway enrichment analysis of the DEmRNAs in the network were performed. In the GO enrichment analysis (Figure 3a, Table S6), the DEmR- ceRNA network were performed. In the GO enrichment analysis (Figure 3a, Table S6), the NAs were involved in the biological processes of enhancer sequence-specific DNA bind- DEmRNAs were involved in the biological processes of enhancer sequence-specific DNA ing, enhancer-binding, protein kinase regulator activity and transcription corepressor ac- binding, enhancer-binding, protein kinase regulator activity and transcription corepressor tivity. The significant KEGG pathway (Figure 3b, Table S7) enrichment included the fork- activity. The significant KEGG pathway (Figure 3b, Table S7) enrichment included the head box O (FoxO) signaling pathway, autophagy and cellular senescence. forkhead box O (FoxO) signaling pathway, autophagy and cellular senescence. Biology 2021, 10, 432 7 of 15 Biology 2021, 10, x 8 of 16 (a) (b) Figure 3. Bar chart of the GO and KEGG analysis of DEmRNAs in the ceRNA network. (a) Bar chart of the significant GO Figure 3. Bar chart of the GO and KEGG analysis of DEmRNAs in the ceRNA network. (a) Bar chart of the significant GO terms with DEmRNAs identified by the clusterProfiler package in R. The color of the bar represents the −log10 transformed terms with DEmRNAs identified by the clusterProfiler package in R. The color of the bar represents the log transformed p values of the GO terms and the abscissa corresponds with the number of genes in the GO terms; (b) bar chart of signifi- p values of the GO terms and the abscissa corresponds with the number of genes in the GO terms; (b) bar chart of significant cant KEGG enriched pathways with DEmRNAs identified by the clusterProfiler package in R. The color of the bar repre- KEGG enriched pathways with DEmRNAs identified by the clusterProfiler package in R. The color of the bar represents sents the −log10 transformed p values of the KEGG enriched pathways and the abscissa corresponds with the number of the log10 transformed p values of the KEGG enriched pathways and the abscissa corresponds with the number of genes in the KEGG enriched pathways. GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; DElncRNAs: differentially expressed lncRNAs; DEmRNAs: differentially expressed mRNAs. Biology 2021, 10, x 9 of 16 Biology 2021, 10, 432 8 of 15 genes in the KEGG enriched pathways. GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; DElncRNAs: differentially expressed lncRNAs; DEmRNAs: differentially expressed mRNAs. 3.4. PPI Network and the Core ceRNA Subnetwork Construction 3.4. PPI Network and the Core ceRNA Subnetwork Construction The combination of the STRING tool and Cytoscape software constructed the PPI The combination of the STRING tool and Cytoscape software constructed the PPI network. A total of 99 nodes and 170 interactions were screened from the PPI network. After network. A total of 99 nodes and 170 interactions were screened from the PPI network. a MCODE analysis, five hub-genes were incorporated into the PPI network (Figure 4a,b). After a MCODE analysis, five hub-genes were incorporated into the PPI network (Figure Based on the ceRNA network, the hub-genes associated with the core ceRNA subnetwork 4a,b). Based on the ceRNA network, the hub-genes associated with the core ceRNA sub- (Tables S8 and S9) were established including 20 miRNA-lncRNA and 7 miRNA-mRNA network (Tables S8 and S9) were established including 20 miRNA-lncRNA and 7 miRNA- edges and 19 nodes (9 lncRNAs (HCG9, POLR2J4, HLA-F-AS1, HCG18, DLEU2, AC074212, mRNA edges and 19 nodes (9 lncRNAs (HCG9, POLR2J4, HLA-F-AS1, HCG18, DLEU2, SNX29P2, LINC00115 and CTC-338M12), 5 miRNAs (hsa-mir-1297, hsa-mir-137, hsa-mir- AC074212, SNX29P2, LINC00115 and CTC-338M12), 5 miRNAs (hsa-mir-1297, hsa-mir- 17-5p, 137, hsa-mir hsa-mir -20b-5p -17-5p, hs and a-mir hsa-mir -20b-5p -25b-5p) and hsand a-mir 5-25b mRNAs -5p) an (AGO4, d 5 mRNA E2F1, s (AG CCND1, O4, E2 TNRC6B F1, CCND1, TNRC6B and ETS1)) (Figure 5). and ETS1)) (Figure 5). (a) (b) Figure 4. PPI network of DEmRNAs and Cluster 1 using the MCODE plug-in from the whole PPI network. (a) PPI network Figure 4. PPI network of DEmRNAs and Cluster 1 using the MCODE plug-in from the whole PPI network. (a) PPI network of DEmRNAs in the ceRNA network obtained from the STRING database. (b) The first cluster is of PPI followed by the of DEmRNAs in the ceRNA network obtained from the STRING database. (b) The first cluster is of PPI followed by the Biology 2021, 10, x 10 of 16 MCODE plug-in of Cytoscape. PPI: Protein-protein interaction; DEmRNAs: differentially expressed mRNAs; MCODE: MCODE plug-in of Cytoscape. PPI: Protein-protein interaction; DEmRNAs: differentially expressed mRNAs; MCODE: molecular complex. molecular complex. Figure 5. The hub-genes associated with the ceRNA network in a cSCC. MRNAs, miRNAs and Figure 5. The hub-genes associated with the ceRNA network in a cSCC. MRNAs, miRNAs and lncRNAs are indicated by the circle in green, rectangle in red and triangle in blue, respectively. lncRNAs are indicated by the circle in green, rectangle in red and triangle in blue, respectively. 3.5. Verification of the DEGs in the Core ceRNA Subnetwork The expression of 9 DElncRNAs and 5 DEmRNAs in the hub-genes associated with the core ceRNA subnetwork were compared between the cSCC and the NS group using the GSE7553 dataset. Finally, one lncRNA HLA-F antisense RNA 1 (HLA-F-AS1) and 3 DEmRNAs (Argonaute-4 (AGO4), Cyclin D1 (CCND1) and E2F transcription factor 1 (E2F1)) were considered to be validated of which the expression patterns were consistent with the results of GSE42677 and GSE45164 (Figure 6). (a) (b) Biology 2021, 10, x 10 of 16 Biology 2021, 10, 432 9 of 15 Figure 5. The hub-genes associated with the ceRNA network in a cSCC. MRNAs, miRNAs and lncRNAs are indicated by the circle in green, rectangle in red and triangle in blue, respectively. 3.5. Verification of the DEGs in the Core ceRNA Subnetwork 3.5. Verification of the DEGs in the Core ceRNA Subnetwork The expression of 9 DElncRNAs and 5 DEmRNAs in the hub-genes associated with The expression of 9 DElncRNAs and 5 DEmRNAs in the hub-genes associated with the core ceRNA subnetwork were compared between the cSCC and the NS group using the core ceRNA subnetwork were compared between the cSCC and the NS group using the GSE7553 dataset. Finally, one lncRNA HLA-F antisense RNA 1 (HLA-F-AS1) and 3 the GSE7553 dataset. Finally, one lncRNA HLA-F antisense RNA 1 (HLA-F-AS1) and 3 DEmRNAsDEm (Argonaute-4 RNAs (Argon (AGO4), auteCyclin -4 (AGO4 D1 (CCND1) ), Cyclin and D1 (CCN E2F transcription D1) and E2F factor transcription 1 factor 1 (E2F1)) were considered to be validated of which the expression patterns were consistent (E2F1)) were considered to be validated of which the expression patterns were consistent with the results of GSE42677 and GSE45164 (Figure 6). with the results of GSE42677 and GSE45164 (Figure 6). Biology 2021, 10, x 11 of 16 (a) (b) (c) (d) Figure 6. External validation of the four core genes in the GSE7553 datasets. (a) The relative Figure 6. External validation of the four core genes in the GSE7553 datasets. (a) The relative expression level of HLA-F- AS1 between expr cSCC ession and level NS sam of HLA-F-AS1 ples; (b) the between relativecSCC expres and sion NSle samples; vel of AGO4 (b) thebetwee relative n expr cSCC ession and level NS samples; of (c) the relative expression level of CCND1 between cSCC and NS samples; (d) the relative expression level of E2F1 between cSCC AGO4 between cSCC and NS samples; (c) the relative expression level of CCND1 between cSCC and NS samples. The log2 transformed expression level corresponded with the relative expression level. Data are presented and NS samples; (d) the relative expression level of E2F1 between cSCC and NS samples. The log as a median with interquartile ranges. An unpaired t-test was used for the comparison of the relative expression level transformed expression level corresponded with the relative expression level. Data are presented as between the two groups. a median with interquartile ranges. An unpaired t-test was used for the comparison of the relative expression level between the two groups. 4. Discussion A cutaneous squamous cell carcinoma (cSCC) derived from the malignant prolifera- tion of keratinocytes is the second most frequent cause of non-melanoma skin cancer [31]. However, a cSCC usually displays a primary cSCC with benign clinical behavior and can be successfully treated with surgery. It may progress to be locally invasive and metastatic [32]. The etiology and pathogenesis of a cSCC have not been well elucidated. Recently, accumulating studies have emphasized the regulatory role of lncRNAs as ceRNAs in the occurrence and progression of many types of cancers including cSCCs [15–18]. We then constructed an lncRNA-associated ceRNA network based on the transcriptome expres- sion profiles from the GEO database to identify new targets with a potential diagnostic or curative value for the cSCC. In total, 24 DElncRNAs and 3221 DEmRNAs were identified between the cSCC tu- mors compared with NS samples. The DElncRNAs and DEmRNAs with a strong relation- ship were incorporated with their common targeted miRNAs to construct the ceRNA net- work. The DEmRNAs in the ceRNA network were then subjected to the GO enrichment and KEGG pathway enrichment analyses. The pathway analysis showed that the AMPK signaling pathway [33] and autophagy [34] participate in the pathogenesis of a cSCC. We established a hub-gene-associated core ceRNA subnetwork including nine lncRNAs, five predicted miRNAs (hsa-mir-1297, hsa-mir-137, hsa-mir-17-5p, hsa-mir-20b- 5p and hsa-mir-25b-5p) and five mRNAs (AGO4, E2F1, CCND1, TNRC6B and ETS1). Three DEmRNAs (E2F1, CCND1 and AGO4) were validated in the independent dataset GSE7553. E2F1 and CCND1 are both cell cycle regulatory proteins. The upregulated ex- pression of E2F1 and its family number E2F7 in cSCCs may contribute to cSCC formation [35]. E2F1 was shown to inhibit the differentiation of keratinocytes at the terminal phase and inactivate MMP [36,37]. The aberrant activation of E2F1 also participated in the cellu- lar proliferation, differentiation and apoptosis of colon cancer and prostate cancer and is often related to a poor prognosis [38,39]. CCND1 has been discovered to be upregulated in SCCs such as laryngeal squamous cell carcinomas (LSCCs), head and neck squamous cell carcinomas (HNSCCs) and oral squamous cell carcinomas (OSCCs) [40–42]. As for the AGO4, it may exert a function as an effector protein with DNA methylation via RNA- Biology 2021, 10, 432 10 of 15 4. Discussion A cutaneous squamous cell carcinoma (cSCC) derived from the malignant prolifera- tion of keratinocytes is the second most frequent cause of non-melanoma skin cancer [31]. However, a cSCC usually displays a primary cSCC with benign clinical behavior and can be successfully treated with surgery. It may progress to be locally invasive and metastatic [32]. The etiology and pathogenesis of a cSCC have not been well elucidated. Recently, accumu- lating studies have emphasized the regulatory role of lncRNAs as ceRNAs in the occurrence and progression of many types of cancers including cSCCs [15–18]. We then constructed an lncRNA-associated ceRNA network based on the transcriptome expression profiles from the GEO database to identify new targets with a potential diagnostic or curative value for the cSCC. In total, 24 DElncRNAs and 3221 DEmRNAs were identified between the cSCC tumors compared with NS samples. The DElncRNAs and DEmRNAs with a strong relationship were incorporated with their common targeted miRNAs to construct the ceRNA network. The DEmRNAs in the ceRNA network were then subjected to the GO enrichment and KEGG pathway enrichment analyses. The pathway analysis showed that the AMPK signaling pathway [33] and autophagy [34] participate in the pathogenesis of a cSCC. We established a hub-gene-associated core ceRNA subnetwork including nine lncR- NAs, five predicted miRNAs (hsa-mir-1297, hsa-mir-137, hsa-mir-17-5p, hsa-mir-20b-5p and hsa-mir-25b-5p) and five mRNAs (AGO4, E2F1, CCND1, TNRC6B and ETS1). Three DEmRNAs (E2F1, CCND1 and AGO4) were validated in the independent dataset GSE7553. E2F1 and CCND1 are both cell cycle regulatory proteins. The upregulated expression of E2F1 and its family number E2F7 in cSCCs may contribute to cSCC formation [35]. E2F1 was shown to inhibit the differentiation of keratinocytes at the terminal phase and inactivate MMP [36,37]. The aberrant activation of E2F1 also participated in the cellular proliferation, differentiation and apoptosis of colon cancer and prostate cancer and is often related to a poor prognosis [38,39]. CCND1 has been discovered to be upregulated in SCCs such as laryngeal squamous cell carcinomas (LSCCs), head and neck squamous cell carcinomas (HNSCCs) and oral squamous cell carcinomas (OSCCs) [40–42]. As for the AGO4, it may exert a function as an effector protein with DNA methylation via RNA-dependent DNA methylation (RdDM) [43]. A DNA methyltransferase 3A (DNMT3A)/AGO4 complex abol- ished miR-181-5p inhibition of gene expression via the cytosine methylation of miR-181-5p and resulted in an aggressive outcome of a glioblastoma multiforme (GBM) [44]. Among the five predicted miRNAs, miR-17-5p had been validated to be highly and significantly upregulated in a cSCC by a real-time polymerase chain reaction (RT-PCR) between cSCC tissues and NS samples [45]. MiR-17-5p and miR-20b-5p are both members of the miR-17-92 family. Moreover, the interaction of miR-17 with CCND1 and E2F1 has been studied. In a study of a lung adenocarcinoma, miR-20b-5p was found to directly and independently inhibit the activity of reporter genes and decrease the CCND1 protein 0 0 expression [46]. MiR-20b-5p could sponge with the 3 untranslated region (3 UTR) of CCND1, repressing the CCND1 expression in a lung cancer cell [46]. MiR-20b-5p also targeted with E2F1, inducing a myoblast differentiation and repressing its proliferation [47]. Only one lncRNA HLA-F-AS1 in the core ceRNA subnetwork was validated in the external datasets. Several studies had revealed the dysregulation of HLA-F-AS1 in can- cers [48–50] while the transcript expression of HLA-F-AS1 was verified to be downregulated in lung adenocarcinoma (LAD) samples [49]. It was over-expressed in colorectal cancer (CRC) tissues and CRC cell lines. Huang et al. found that HLA-F-AS1 may participate in the progression of CRC via the miR-330-3p/profilin 1 (PFN1) axis in CRC cells [48]. The elevated HLA-F-AS1 expression induced by STAT3 promoted cell growth and stemness characteristics via sponging with miR-541-3p to upregulate the expression of a TraB domain (TRABD) containing triple-negative breast cancer cells [50]. In summary, according to the predicted core ceRNA subnetwork, lncRNA HLA-F-AS1 may sponge miR-17-5p or miR-20b-5p to regulate the expression of CCND1 and E2F1 in the cSCC. The expression of HLA-F-AS1 in the cSCC samples and its role mechanism in Biology 2021, 10, 432 11 of 15 the biological process of a cSCC as a ceRNA required further research. To highlight the main findings of the validated lncRNAs and miRNAs in the core ceRNA subnetwork, we followed the model presented by Garofoli et al. [51] and made a table where a systematic review of the roles and mechanisms of HLA-F-AS1, miR-17-5p and miR-20b-5p in the pathogenesis of cancers was carried out (Table 1). Table 1. Summary of the roles and mechanisms of discussed lncRNAs and miRNAs in the pathogenesis of cancers. LncRNA/miRNA Gene Targets Mechanisms in Cancers Promotion of cell proliferation and stemness in vitro and tumor growth in vivo by MiR-541-3p the miR-541-3p/TRABD axis in triple-negative breast cancer [50]. MiR-330-3p HLA-F-AS1 Promotion of cell proliferation, migration and invasion and repression of apoptosis in vitro via the miR-330-3p/ PFN1 axis in colorectal cancer [48]. 0 0 3 UTR of CCND1 Repression of the CCND1 protein expression via sponging with the CCND1 3 -UTR E2F1 region in a lung cancer cell [46]. MiR-17-5p Promotion of cell growth, proliferation, migration and colony formation in vitro in a non-metastatic hepatocellular carcinoma by regulating E2F1 expression [52]. Promotion of myoblast differentiation and repression of myoblast proliferation by 3 UTR of CCND1 regulating E2F1 expression [47]. 3 UTR of E2F1 Promotion of cell proliferation, G1/S transition and DNA synthesis by downregulating CCND1 and E2F1 expression in multipotent stromal cells [53]. MiR-20b-5p Promotion of cell migration and proliferation and repression of apoptosis of breast cancer stem cells and promoting tumor growth in vivo by bidirectionally regulating the protein levels of CCND1 and E2F1 in breast cancer [54]. Inhibition of the cell cycle, migration and invasion in vitro and the tumorigenesis in vivo by negatively regulating CCND1 in colon cancer [55]. HLA-F-AS1: HLA-F antisense RNA 1; CCND1: Cyclin D1; E2F1: E2F transcription factor 1; TRABD: TraB domain containing; PFN1: 0 0 profilin 1; 3 UTR: 3 untranslated regions. 5. Conclusions In this study, the cSCC-related lncRNA-miRNA-mRNA ceRNA network was es- tablished and one lncRNA (HLA-F-AS1) and three mRNAs (AGO4, CCND1 and E2F1) in the core ceRNA network were validated. The HLA-F-AS1-miR-17-5p/miR-20b-5p- CCND1/E2F1 axis may take part in the development of a cSCC. However, there are several limitations in our research that should be considered. First, miRNA expression datasets are needed to verify the network. Second, a relatively large sample size should be collected to assure the results are more reliable in future studies. Moreover, the relation between the ceRNA network and clinical features could not be established due to the absence of related clinical information. Our research provides a new therapeutic idea and regulatory roles of ceRNAs in the pathogenesis of a cSCC. They may be candidate targets with a promising diagnostic or therapeutic value with a cSCC. Further functional investigations mediated by the ceRNA would have to be verified. Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/biology10050432/s1, Figure S1: Data preprocessing. Table S1: The lncRNA-mRNA pairs with significant correlations, Table S2: DElncRNA-miRNA pairs of the ceRNA network, Table S3: DEmRNA-miRNA pairs of the ceRNA network, Table S4: Over-expressed and down-expressed lncR- NAs of the ceRNA network, Table S5: Over-expressed and down-expressed mRNAs of the ceRNA network, Table S6: GO enrichment terms of DEmRNAs in the ceRNA network, Table S7: KEGG pathway enrichment of DEmRNAs in the ceRNA network, Table S8: Details of DElncRNA-miRNA pairs of the core ceRNA network, Table S9: DEmRNA-miRNA pairs of the core ceRNA network. Author Contributions: Y.X., Y.D. (Yunhua Deng) and Q.Q. conceived the idea. Y.X., Y.D. (Yingying Dong), M.W., H.L., Q.S., Q.G. and Q.Q. performed the literature review, data analysis and wrote the manuscript. Y.D. (Yunhua Deng) and Q.Q. revised the manuscript and helped to improve the quality of the manuscript. All authors have read and agreed to the published version of the manuscript. Biology 2021, 10, 432 12 of 15 Funding: This research was funded by the National Natural Science Foundation of China (Grant No. 81773306). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Data were obtained from the GEO database and are available at: www.ncbi.nlm.nih.gov/geo (accessed on 24 August 2020) with the permission of the GEO database. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations SCC Squamous Cell Carcinoma cSCC Cutaneous Squamous Cell Carcinoma SCCIS Squamous Cell Carcinoma (SCC) In Situ NcRNAs Non-Coding RNAs LncRNAs Long Non-Coding RNAs MiRNAs MicroRNAs CeRNAs Competing Endogenous RNAs NS Normal Skin DEGs Differentially Expressed Genes DElncRNAs Differentially Expressed LncRNAs DEmRNAs Differentially Expressed MRNAs GEO Gene Expression Omnibus GO Gene Ontology BP Biological Processes KEGG Kyoto Encyclopedia of Genes and Genomes STRING Search Tool for the Retrieval of Interacting Genes MCODE Molecular Complex Detection FDR False Discovery Rate FC Fold Change PCA Principal Component Analysis PCC Pearson’s Correlation Coefficient CDC20 Cell Division Cycle 20 EGFR Epidermal Growth Factor Receptor PRAF2 Prenylated Rab Acceptor 1 Domain Family Member 2 FoxO Forkhead box O HLA-F-AS1 HLA-F Antisense RNA 1 CCND1 Cyclin D1 E2F1 E2F Transcription Factor 1 AGO4 Argonaute-4 LAD Lung Adenocarcinoma CRC Colorectal Cancer PFN1 Profilin 1 ESCC Esophageal Squamous Cell Carcinoma TSCC Tongue Squamous Cell Carcinoma Progression LSCC Laryngeal Squamous Cell Carcinoma HNSCC Head and Neck Squamous Cell Carcinoma OSCC Oral Squamous Cell Carcinoma RdDM RNA-Dependent DNA Methylation DNMT3A DNA Methyltransferase 3A GBM Glioblastoma Multiforme AK Actinic Keratosis NEAT1 Nuclear Paraspeckle Assembly Transcript 1 RT-PCR Real-time Polymerase Chain Reaction 0 0 3 UTRs 3 Untranslated Regions TRABD TraB Domain Containing Biology 2021, 10, 432 13 of 15 References 1. Motaparthi, K.; Kapil, J.P.; Velazquez, E.F. Cutaneous Squamous Cell Carcinoma: Review of the Eighth Edition of the American Joint Committee on Cancer Staging Guidelines, Prognostic Factors, and Histopathologic Variants. Adv. Anat. Pathol. 2017, 24, 171–194. [CrossRef] 2. Waldman, A.; Schmults, C. Cutaneous Squamous Cell Carcinoma. Hematol. Oncol. Clin. N. Am. 2019, 33, 1–12. [CrossRef] 3. Fernandez Figueras, M.T. From actinic keratosis to squamous cell carcinoma: Pathophysiology revisited. J. Eur. Acad. Derm. Venereol. 2017, 31, 5–7. [CrossRef] [PubMed] 4. Soura, E.; Gagari, E.; Stratigos, A. Advanced cutaneous squamous cell carcinoma: How is it defined and what new therapeutic approaches are available. Curr. Opin. Oncol. 2019, 31, 461–468. [CrossRef] [PubMed] 5. Esteller, M. Non-coding RNAs in human disease. Nat. Rev. Genet. 2011, 12, 861–874. [CrossRef] 6. Bai, X.; Zhou, Y.; Chen, P.; Yang, M.; Xu, J. MicroRNA-142-5p induces cancer stem cell-like properties of cutaneous squamous cell carcinoma via inhibiting PTEN. J. Cell. Biochem. 2018, 119, 2179–2188. [CrossRef] [PubMed] 7. Wang, Y.; Deng, X.; Dai, Y.; Niu, X.; Zhou, M. miR-27a Downregulation Promotes Cutaneous Squamous Cell Carcinoma Progression via Targeting EGFR. Front. Oncol. 2020, 9, 1565. [CrossRef] 8. Quinn, J.J.; Chang, H.Y. Unique features of long non-coding RNA biogenesis and function. Nat. Rev. Genet. 2016, 17, 47–62. [CrossRef] 9. Rotondo, J.C.; Selvatici, R.; Di Domenico, M.; Marci, R.; Vesce, F.; Tognon, M.; Martini, F. Methylation loss at H19 imprinted gene correlates with methylenetetrahydrofolate reductase gene promoter hypermethylation in semen samples from infertile males. Epigenetics 2013, 8, 990–997. [CrossRef] 10. Di Martino, M.T.; Riillo, C.; Scionti, F.; Grillone, K.; Polerà, N.; Caracciolo, D.; Arbitrio, M.; Tagliaferri, P.; Tassone, P. miRNAs and lncRNAs as Novel Therapeutic Targets to Improve Cancer Immunotherapy. Cancers (Basel) 2021, 13, 1587. [CrossRef] 11. Mei, X.L.; Zhong, S. Long noncoding RNA LINC00520 prevents the progression of cutaneous squamous cell carcinoma through the inactivation of the PI3K/Akt signaling pathway by downregulating EGFR. Chin. Med. J. (Engl.) 2019, 132, 454–465. [CrossRef] 12. Piipponen, M.; Heino, J.; Kähäri, V.M.; Nissinen, L. Long non-coding RNA PICSAR decreases adhesion and promotes migration of squamous carcinoma cells by downregulating 2 1 and 5 1 integrin expression. Biol. Open 2018, 7, bio037044. [CrossRef] 13. Yamamura, S.; Imai-Sumida, M.; Tanaka, Y.; Dahiya, R. Interaction and cross-talk between non-coding RNAs. Cell. Mol. Life Sci. 2018, 75, 467–484. [CrossRef] 14. López-Urrutia, E.; Bustamante Montes, L.P.; Ladrón de Guevara Cervantes, D.; Pérez-Plasencia, C.; Campos-Parra, A.D. Crosstalk Between Long Non-coding RNAs, Micro-RNAs and mRNAs: Deciphering Molecular Mechanisms of Master Regulators in Cancer. Front. Oncol. 2019, 9, 669. [CrossRef] 15. Li, Y.; Chen, D.; Gao, X.; Li, X.; Shi, G. LncRNA NEAT1 Regulates Cell Viability and Invasion in Esophageal Squamous Cell Carcinoma through the miR-129/CTBP2 Axis. Dis. Markers 2017, 9, 5314649. [CrossRef] 16. Ding, J.; Yang, C.; Yang, S. LINC00511 interacts with miR-765 and modulates tongue squamous cell carcinoma progression by targeting LAMC2. J. Oral. Pathol. Med. 2018, 47, 468–476. [CrossRef] [PubMed] 17. Zhao, J.; Bai, X.; Feng, C.; Shang, X.; Xi, Y. Long Non-Coding RNA HCP5 Facilitates Cell Invasion And Epithelial-Mesenchymal Transition In Oral Squamous Cell Carcinoma By miR-140-5p/SOX4 Axis. Cancer Manag. Res. 2019, 11, 10455–10462. [CrossRef] [PubMed] 18. Yu, G.J.; Sun, Y.; Zhang, D.W.; Zhang, P. Long non-coding RNA HOTAIR functions as a competitive endogenous RNA to regulate PRAF2 expression by sponging miR-326 in cutaneous squamous cell carcinoma. Cancer Cell Int. 2019, 19, 270. [CrossRef] 19. Johnson, W.E.; Li, C.; Rabinovic, A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 2007, 8, 118–127. [CrossRef] [PubMed] 20. Konishi, T. Principal component analysis for designed experiments. BMC Bioinform. 2015, 16, S7. [CrossRef] 21. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015, 43, e47. [CrossRef] [PubMed] 22. Wickham, H. ggplot2: Elegant Graphics for Data Analysis; Springer: New York, NY, USA, 2016. 23. Jeggari, A.; Marks, D.S.; Larsson, E. miRcode: A map of putative microRNA target sites in the long non-coding transcriptome. Bioinformatics 2012, 28, 2062–2063. [CrossRef] 24. Chou, C.H.; Shrestha, S.; Yang, C.D.; Chang, N.W.; Lin, Y.L.; Liao, K.W.; Huang, W.C.; Sun, T.H.; Tu, S.J.; Lee, W.H.; et al. miRTarBase update 2018: A resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 2018, 46, D296–D302. [CrossRef] [PubMed] 25. Chen, Y.; Wang, X. miRDB: An online database for prediction of functional microRNA targets. Nucleic Acids Res. 2020, 48, D127–D131. [CrossRef] [PubMed] 26. Lewis, B.P.; Shih, I.-H.; Jones-Rhoades, M.W.; Bartel, D.P.; Burge, C.B. Prediction of mammalian microRNA targets. Cell 2003, 115, 787–798. [CrossRef] 27. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [CrossRef] [PubMed] 28. Yu, G.; Wang, L.-G.; Han, Y.; He, Q.-Y. clusterProfiler: An R package for comparing biological themes among gene clusters. Omics 2012, 16, 284–287. [CrossRef] Biology 2021, 10, 432 14 of 15 29. Szklarczyk, D.; Franceschini, A.; Wyder, S.; Forslund, K.; Heller, D.; Huerta-Cepas, J.; Simonovic, M.; Roth, A.; Santos, A.; Tsafou, K.P.; et al. STRING v10: Protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015, 43, D447–D452. [CrossRef] 30. Bader, G.D.; Hogue, C.W.V. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinform. 2003, 4, 2. [CrossRef] 31. Que, S.K.T.; Zwald, F.O.; Schmults, C.D. Cutaneous squamous cell carcinoma: Management of advanced and high-stage tumors. J. Am. Acad. Derm. 2018, 78, 249–261. [CrossRef] 32. Brougham, N.D.; Dennett, E.R.; Cameron, R.; Tan, S.T. The incidence of metastasis from cutaneous squamous cell carcinoma and the impact of its risk factors. J. Surg. Oncol. 2012, 106, 811–815. [CrossRef] 33. Piipponen, M.; Nissinen, L.; Farshchian, M.; Riihilä, P.; Kivisaari, A.; Kallajoki, M.; Peltonen, J.; Peltonen, S.; Kähäri, V.M. Long Noncoding RNA PICSAR Promotes Growth of Cutaneous Squamous Cell Carcinoma by Regulating ERK1/2 Activity. J. Investig. Derm. 2016, 136, 1701–1710. [CrossRef] [PubMed] 34. Zhou, W.; Zhang, S.; Li, J.; Li, Z.; Wang, Y.; Li, X. lncRNA TINCR participates in ALA-PDT-induced apoptosis and autophagy in cutaneous squamous cell carcinoma. J. Cell. Biochem. 2019, 120, 13893–13902. [CrossRef] 35. Endo-Munoz, L.; Dahler, A.; Teakle, N.; Rickwood, D.; Hazar-Rethinam, M.; Abdul-Jabbar, I.; Sommerville, S.; Dickinson, I.; Kaur, P.; Paquet-Fifield, S.; et al. E2F7 can regulate proliferation, differentiation, and apoptotic responses in human keratinocytes: Implications for cutaneous squamous cell carcinoma formation. Cancer Res. 2009, 69, 1800–1808. [CrossRef] [PubMed] 36. Das Mahapatra, K.; Pasquali, L.; Søndergaard, J.N.; Lapins, J.; Nemeth, I.B.; Baltás, E.; Kemény, L.; Homey, B.; Moldovan, L.-I.; Kjems, J.; et al. A comprehensive analysis of coding and non-coding transcriptomic changes in cutaneous squamous cell carcinoma. Sci. Rep. 2020, 10, 3637. [CrossRef] [PubMed] 37. Nagarajan, P.; Chin, S.S.; Wang, D.; Liu, S.; Sinha, S.; Garrett-Sinha, L.A. Ets1 blocks terminal differentiation of keratinocytes and induces expression of matrix metalloproteases and innate immune mediators. J. Cell Sci. 2010, 123, 3566–3575. [CrossRef] 38. Chun, J.N.; Cho, M.; Park, S.; So, I.; Jeon, J.-H. The conflicting role of E2F1 in prostate cancer: A matter of cell context or interpretational flexibility? Biochim. Biophys. Acta Rev. Cancer 2020, 1873, 188336. [CrossRef] 39. Fang, Z.; Lin, M.; Li, C.; Liu, H.; Gong, C. A comprehensive review of the roles of E2F1 in colon cancer. Am. J. Cancer Res. 2020, 10, 757–768. 40. Hanken, H.; Gröbe, A.; Cachovan, G.; Smeets, R.; Simon, R.; Sauter, G.; Heiland, M.; Blessmann, M. CCND1 amplification and cyclin D1 immunohistochemical expression in head and neck squamous cell carcinomas. Clin. Oral Investig. 2014, 18, 269–276. [CrossRef] 41. Wang, X.; Wang, C.; Yan, G.; Kang, Y.; Sun, G.; Wang, S.; Zou, R.; Sun, H.; Zeng, K.; Song, H.; et al. BAP18 is involved in upregulation of CCND1/2 transcription to promote cell growth in oral squamous cell carcinoma. EBioMedicine 2020, 53, 102685. [CrossRef] 42. Zang, Y.; Li, J.; Wan, B.; Tai, Y. circRNA circ-CCND1 promotes the proliferation of laryngeal squamous cell carcinoma through elevating CCND1 expression via interacting with HuR and miR-646. J. Cell. Mol. Med. 2020, 24, 2423–2433. [CrossRef] 43. Chalertpet, K.; Pin-On, P.; Aporntewan, C.; Patchsung, M.; Ingrungruanglert, P.; Israsena, N.; Mutirangura, A. Argonaute 4 as an Effector Protein in RNA-Directed DNA Methylation in Human Cells. Front. Genet. 2019, 10, 645. [CrossRef] 44. Cheray, M.; Etcheverry, A.; Jacques, C.; Pacaud, R.; Bougras-Cartron, G.; Aubry, M.; Denoual, F.; Peterlongo, P.; Nadaradjane, A.; Briand, J.; et al. Cytosine methylation of mature microRNAs inhibits their functions and is associated with poor prognosis in glioblastoma multiforme. Mol. Cancer 2020, 19, 36. [CrossRef] 45. Sand, M.; Hessam, S.; Amur, S.; Skrygan, M.; Bromba, M.; Stockfleth, E.; Gambichler, T.; Bechara, F.G. Expression of oncogenic miR-17-92 and tumor suppressive miR-143-145 clusters in basal cell carcinoma and cutaneous squamous cell carcinoma. J. Derm. Sci. 2017, 86, 142–148. [CrossRef] 46. Ding, H.; Luo, Y.; Hu, K.; Liu, P.; Xiong, M. Linc00467 promotes lung adenocarcinoma proliferation via sponging miR-20b-5p to activate CCND1 expression. Onco Targets Ther. 2019, 12, 6733–6743. [CrossRef] 47. Luo, W.; Li, G.; Yi, Z.; Nie, Q.; Zhang, X. E2F1-miR-20a-5p/20b-5p auto-regulatory feedback loop involved in myoblast proliferation and differentiation. Sci. Rep. 2016, 6, 27904. [CrossRef] [PubMed] 48. Huang, Y.; Sun, H.; Ma, X.; Zeng, Y.; Pan, Y.; Yu, D.; Liu, Z.; Xiang, Y. HLA-F-AS1/miR-330-3p/PFN1 axis promotes colorectal cancer progression. Life Sci. 2020, 254, 117180. [CrossRef] 49. Yang, Z.; Li, H.; Wang, Z.; Yang, Y.; Niu, J.; Liu, Y.; Sun, Z.; Yin, C. Microarray expression profile of long non-coding RNAs in human lung adenocarcinoma. Thorac. Cancer 2018, 9, 1312–1322. [CrossRef] [PubMed] 50. Wu, D.; Jia, H.; Zhang, Z.; Li, S. STAT3-induced HLA-F-AS1 promotes cell proliferation and stemness characteristics in triple negative breast cancer cells by upregulating TRABD. Bioorg. Chem. 2021, 109, 104722. [CrossRef] [PubMed] 51. Garofoli, M.; Volpicella, M.; Guida, M.; Porcelli, L.; Azzariti, A. The Role of Non-Coding RNAs as Prognostic Factor, Predictor of Drug Response or Resistance and Pharmacological Targets, in the Cutaneous Squamous Cell Carcinoma. Cancers (Basel) 2020, 12, 2552. [CrossRef] 52. El Tayebi, H.M.; Omar, K.; Hegy, S.; El Maghrabi, M.; El Brolosy, M.; Hosny, K.A.; Esmat, G.; Abdelaziz, A.I. Repression of miR-17-5p with elevated expression of E2F-1 and c-MYC in non-metastatic hepatocellular carcinoma and enhancement of cell growth upon reversing this expression pattern. Biochem. Biophys. Res. Commun. 2013, 434, 421–427. [CrossRef] [PubMed] Biology 2021, 10, 432 15 of 15 53. Tai, L.; Huang, C.J.; Choo, K.B.; Cheong, S.K.; Kamarul, T. Oxidative Stress Down-Regulates MiR-20b-5p, MiR-106a-5p and E2F1 Expression to Suppress the G1/S Transition of the Cell Cycle in Multipotent Stromal Cells. Int. J. Med. Sci. 2020, 17, 457–470. [CrossRef] [PubMed] 54. Xia, L.; Li, F.; Qiu, J.; Feng, Z.; Xu, Z.; Chen, Z.; Sun, J. Oncogenic miR-20b-5p contributes to malignant behaviors of breast cancer stem cells by bidirectionally regulating CCND1 and E2F1. BMC Cancer 2020, 20, 949. [CrossRef] 55. Yang, H.; Lin, J.; Jiang, J.; Ji, J.; Wang, C.; Zhang, J. miR-20b-5p functions as tumor suppressor microRNA by targeting cyclinD1 in colon cancer. Cell Cycle 2020, 19, 2939–2954. [CrossRef] [PubMed] http://www.deepdyve.com/assets/images/DeepDyve-Logo-lg.png Biology Multidisciplinary Digital Publishing Institute

Identifying an lncRNA-Related ceRNA Network to Reveal Novel Targets for a Cutaneous Squamous Cell Carcinoma

Loading next page...
 
/lp/multidisciplinary-digital-publishing-institute/identifying-an-lncrna-related-cerna-network-to-reveal-novel-targets-OeADj9Z964

References (57)

Publisher
Multidisciplinary Digital Publishing Institute
Copyright
© 1996-2021 MDPI (Basel, Switzerland) unless otherwise stated Disclaimer The statements, opinions and data contained in the journals are solely those of the individual authors and contributors and not of the publisher and the editor(s). MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Terms and Conditions Privacy Policy
ISSN
2079-7737
DOI
10.3390/biology10050432
Publisher site
See Article on Publisher Site

Abstract

biology Article Identifying an lncRNA-Related ceRNA Network to Reveal Novel Targets for a Cutaneous Squamous Cell Carcinoma 1 1 1 , 2 , 3 3 1 Yaqin Xu , Yingying Dong , Yunhua Deng *, Qianrong Qi * , Mi Wu , Hongmei Liang , Qiuyun She and Qing Guo Department of Dermatology, Tongji Medical College, Huazhong University of Science and Technology, Wuhan 430030, China; xuyaqin@hust.edu.cn (Y.X.); yingyingdong@hust.edu.cn (Y.D.); qiuyunshe@hust.edu.cn (Q.S.) Department of Obstetrics & Gynecology, University of California, Irvine, CA 92697, USA Department of Immunology, School of Basic Medicine, Tongji Medical College, Huazhong University of Science and Technology, Wuhan 430030, China; wumi@hust.edu.cn (M.W.); lianghongmei@hust.edu.cn (H.L.); qingguo@hust.edu.cn (Q.G.) * Correspondence: yhdeng@tjh.tjmu.edu.cn (Y.D.); qianronq@hs.uci.edu (Q.Q.) Simple Summary: The exact functions and molecular mechanism of lncRNAs, acting as competitive endogenous RNAs in a cutaneous squamous cell carcinoma, remain unexplored. The present study was conducted to identify the differentially expressed lncRNAs and mRNAs and establish the lncRNA-related competing endogenous RNA networks associated with a cutaneous squamous cell carcinoma. A competing endogenous RNA network consisting of 137 miRNA-lncRNA and 221 miRNA-mRNA pairs was constructed. As for the functional analysis of the mRNAs in the network, a FoxO signaling pathway, an autophagy and cellular senescence were the top enrichment terms based on the Kyoto Encyclopedia of Genes and Genomes analysis. We identified five core Citation: Xu, Y.; Dong, Y.; Deng, Y.; mRNAs and built a core mRNA-associated competing endogenous RNA network. Finally, one Qi, Q.; Wu, M.; Liang, H.; She, Q.; lncRNA HLA-F-AS1 and three mRNAs named AGO4, E2F1 and CCND1 in the core mRNA-associated Guo, Q. Identifying an lncRNA- competing endogenous RNA network were validated with the same expression patterns. The core Related ceRNA Network to Reveal mRNAs and their associated lncRNAs may provide potential therapeutic targets for cutaneous Novel Targets for a Cutaneous squamous cell carcinomas. Squamous Cell Carcinoma. Biology 2021, 10, 432. https://doi.org/ Abstract: A cutaneous squamous cell carcinoma (cSCC) derived from keratinocytes is the second 10.3390/biology10050432 most common cause of non-melanoma skin cancer. The accumulation of the mutational burden of genes and cellular DNA damage caused by the risk factors (e.g., exposure to ultraviolet radiation) Academic Editor: Loredana Moro contribute to the aberrant proliferation of keratinocytes and the formation of a cSCC. A cSCC Received: 22 March 2021 encompasses a spectrum of diseases that range from recursor actinic keratosis (AK) and squamous Accepted: 9 May 2021 cell carcinoma (SCC) in situ (SCCIS) to invasive cSCCs and further metastatic SCCs. Emerging Published: 13 May 2021 evidence has revealed that lncRNAs are involved in the biological process of a cSCC. According to the ceRNA regulatory theory, lncRNAs act as natural miRNA sponges and interact with miRNA Publisher’s Note: MDPI stays neutral response elements, thereby regulating the mRNA expression of their down-stream targets. This with regard to jurisdictional claims in study was designed to search for the potential lncRNAs that may become potential therapeutic published maps and institutional affil- targets or biomarkers of a cSCC. Considering the spirit of the study to be adequately justified, we iations. collected microarray-based datasets of 19 cSCC tissues and 12 normal skin samples from the GEO database (GSE42677 and GSE45164). After screening the differentially expressed genes via a limma package, we identified 24 differentially expressed lncRNAs (DElncRNAs) and 3221 differentially expressed mRNAs (DEmRNAs). The miRcode, miRTarBase, miRDB and TargetScan databases were Copyright: © 2021 by the authors. used to predict miRNAs that could interact with DElncRNAs and DEmRNAs. A total of 137 miRNA- Licensee MDPI, Basel, Switzerland. lncRNA and 221 miRNA-mRNA pairs were retained in the ceRNA network, consisting of 31 miRNAs, This article is an open access article 11 DElncRNAs and 155 DEmRNAs. For the functional analysis, the top enriched biological process distributed under the terms and was enhancer sequence-specific DNA binding in Gene Ontology (GO) terms. The FoxO signaling conditions of the Creative Commons pathway, autophagy and cellular senescence were the top enrichment terms based on a Kyoto Attribution (CC BY) license (https:// Encyclopedia of Genes and Genomes (KEGG) analysis. The combination of a STRING tool and creativecommons.org/licenses/by/ Cytoscape software (plug-in MCODE) identified five core mRNAs and built a core mRNA-associated 4.0/). Biology 2021, 10, 432. https://doi.org/10.3390/biology10050432 https://www.mdpi.com/journal/biology Biology 2021, 10, 432 2 of 15 ceRNA network. The expression for five identified core mRNAs and their related nine lncRNAs was validated using the external dataset GSE7553. Finally, one lncRNA HLA-F-AS1 and three mRNAs named AGO4, E2F1 and CCND1 were validated with the same expression patterns. We speculate that lncRNA HLA-F-AS1 may sponge miR-17-5p or miR-20b-5p to regulate the expression of CCND1 and E2F1 in the cSCC. The present study may provide potential diagnostic and therapeutic targets for cSCC patients. Keywords: competing endogenous RNAs network; long non-coding RNAs; microRNAs; cutaneous squamous cell carcinoma 1. Introduction Cutaneous squamous cell carcinoma (cSCC), which originates from epidermal ker- atinocytes, ranks as the second most occurring skin cancer [1]. The occurrence of cSCC is due to many risk factors. Beside the most important risk factor, cumulative ultraviolet exposure, other risk factors include old age, immunosuppression, smoking and genetic factors [2]. Moreover, a cSCC commonly has high mutational burden (e.g., TP53, CDKN2A, Ras and NOTCH1). The accumulation of these mutations and cellular DNA damage caused by the risk factors contribute to the uncontrolled proliferation of keratinocytes and the formation of a cSCC and disease progression [3]. Consequently, a cSCC is a spectrum of diseases that commonly arises from precursor actinic keratosis (AK) and a squamous cell carcinoma (SCC) in situ (SCCIS) to an invasive cSCC to a metastatic SCC with a distant spread [4]. Based on the tumor stage, the optimal treatment approach for a cSCC is selected. The cSCC at the early stage is generally identified and can be cured via surgical resection and/or radiotherapy. For cases at the advanced stage (invasive and metastatic cSCCs), a single or a combined therapy, including chemotherapy, radiation therapy and epidermal growth factor receptor (EGFR) inhibitors, have been implemented with limited and unsat- isfactory efficacy with a long-term survival rate of less than 20% [4]. Thus, it is urgent to identify new biomarkers and targets that are more effective for the early diagnosis and treatment of a cSCC. Non-coding RNAs (ncRNAs) incapable of protein coding can be classified according to their length including long non-coding RNAs (lncRNAs) (>200 nucleotides) and small ncRNAs (<200 nucleotides) such as microRNAs (miRNAs), small nucleolar RNAs and transfer RNAs [5]. The dysregulation of the expression of miRNAs has also been elucidated to be implicated in the tumorigenesis of a cSCC. The upregulated miR-142-5p induced cancer stem cell-like properties in a cSCC via the direct regulation of PTEN [6]. MiR-27a directly bound to EGFR, leading to the inactivation of the NF-B pathway and the inhibition of cSCC cell growth [7]. LncRNAs are a wide and diverse class of single-stranded ncRNAs longer than 200 nu- cleotides lacking a protein-coding capacity [8]. LncRNA dysregulation has also been linked to the process of cell differentiation and growth and the pathogenesis of many types of cancers [9,10] including cSCCs [11,12]. LINC00520 targeted EFGR to inactivate the PI3K/Akt pathway and restrain the cSCC progress [11]. The lncRNA PICSAR dys- regulated the collagen and fibronectin expression, resulting in the dysregulation of the adhesion and migration of cSCC cells [12]. Emerging evidence demonstrates that activities of lncRNAs and miRNAs are related to each other in different complex mechanisms [13]. One of these mechanisms is the role of lncRNAs as competing endogenous RNAs (ceR- NAs). Specifically, lncRNAs can compete with mRNA to bind miRNA binding sites, thereby negatively regulating miRNAs and their target genes. In other words, lncRNAs act as a sponge for miRNAs abolishing their inhibitory action to target mRNAs [14]. The lncRNA-miRNA-mRNA regulatory network has been verified to be associated with the progression of several types of squamous cell carcinomas (SCCs) [15–17]. LncRNA nuclear paraspeckle assembly transcript 1 (NEAT1) sponged to miR-129, leading to the depression Biology 2021, 10, 432 3 of 15 of CTBP2 and regulated the esophageal squamous cell carcinoma (ESCC) cell viability and invasion [15]. LINC00511 induced the expression of LAMC2 via sponging to miR-765 and modulated a tongue squamous cell carcinoma progression (TSCC) [16]. Zhao et al. indicated that lncRNA HCP5 functioned as a molecular sponge to absorb miR-140-5p in an oral squamous cell carcinoma (OSCC) and upregulated the SOX4 gene expression [17]. As for cSCCs, lncRNA HOTAIR promoted the progression of a cSCC, exerting a function as a ceRNA through the miR-326/Prenylated Rab acceptor 1 domain family member 2 (PRAF2) axis [18]. In this paper, we compared total RNA expression data based on an array between cSCC and normal skin (NS) samples from the National Center of Biotechnology Information (NCBI) Gene Expression Omnibus (GEO) database. Based on the differentially expressed genes (DEGs), we performed a series of down-stream analyses including the construction of an lncRNA-related ceRNA network, a functional enrichment analysis and a protein- protein interaction analysis. We theorize that lncRNA HLA-F-AS1 may act as the sponge of miR-17-5p and miR-20b-5p to regulate the expression of CCND1 and E2F1 in the cSCC. The expression of HLA-F-AS1 in the cSCC samples and its role mechanism in the biological process of a cSCC as a ceRNA required further research. They may be candidate targets with a promising diagnostic or therapeutic value with a cSCC. The remainder of this paper is organized as follows. Section 2 introduces the dataset collection and presents our designed data processing workflow and describes our down- stream analyses and validation methods. The experiment results are described in Section 3. The Discussion based on the experiment results is analyzed in Section 4. Conclusions are drawn in Section 5. 2. Materials and Methods 2.1. Dataset Collection from the GEO Database The expression profiles of lncRNAs and mRNAs between a cSCC and normal skin (NS) from healthy samples were downloaded from the GEO database (www.ncbi.nlm.nih.gov/ geo, accessed on 24 August 2020). The criterion for selecting the cSCC samples was that the clinical and histological diagnosis of biopsies obtained during surgery was a primary cSCC (SCCIS or invasive cSCC). The datasets were recruited with the standard based on the same microarray platform and containing more than three samples in each group. Finally, GSE42677 and GSE45164 dataset bases on the GPL571 platform were obtained including 20 cSCC and 13 NS samples. 2.2. Data Preprocessing and Screening for DElncRNAs and DEmRNAs Based on the annotation information from Affymetrix Human Genome U133A 2.0 Array platform, the probe sequences were quantified to the expression of lncRNAs and mRNAs after mapping to the human GRCh38 reference genome (http://useast.ensembl. org/info/data/ftp/index.html, accessed on 29 August 2020). We integrated all samples from the two datasets (GSE42677 and GSE45164) and applied ComBat [19] to remove the batch-effects of expression values (Figure S1a,b). A principal component analysis (PCA) [20] was applied to assess the independence of the available genes between 20 cSCC and 13 NS samples. According to the PCA result (Figure S1), we removed one cSCC and one control sample affecting the classification. Data of 19 cSCC patients and 12 NS were selected for a subsequent analysis (Figure 1a). Biology 2021, 10, x 4 of 16 Volcano plots were performed for the DEGs using the ggplot2 package [22] (version 3.3.2; Biology 2021, 10, 432 4 of 15 https://cran.r-project.org/web/packages/ggplot2/index.html, accessed on 29 August 2020) in the R software. Pheatmap package (version: 1.0.12; https://cran.r-project.org/web/pack- ages/pheatmap/, accessed 29 August 2020) in R was used to generate the hierarchical cluster heatmap representing the expression intensity and direction (Figure 1b). (a) Biology 2021, 10, x 5 of 16 (b) (c) Figure 1. Differentially expressed analysis between cSCC and NS tissues from GSE42677 and GSE45164 datasets. (a) The principal component analysis (PCA) of the selected 19 cSCC and 12 NS tissues; (b) volcano plots showing DEGs. The red Figure 1. Differentially expressed analysis between cSCC and NS tissues from GSE42677 and GSE45164 datasets. (a) The dots represent an upregulated expression and the blue dots represent downregulated DEGs; (c) cluster analysis (heatmap) between 19 cSCC and 12 NS tissues. The red stripes display upregulation while the blue stripes display the downregula- principal component analysis (PCA) of the selected 19 cSCC and 12 NS tissues; (b) volcano plots showing DEGs. The red tion of genes. DEGs: differential expression genes. dots represent an upregulated expression and the blue dots represent downregulated DEGs; (c) cluster analysis (heatmap) 2.3. Construction of the cSCC-Associated lncRNA-miRNA-mRNA Network between 19 cSCC and 12 NS tissues. The red stripes Due to displ the absence ay upr of egulation miRNA information while , the the miRco blue de data stripes base [23] display the downregulation (http://www.mircode.org/, accessed on 1 September 2020) was used to predict the targeted miRNAs of DElncRNAs. According to the miRTarBase [24] (http://mirtar- of genes. DEGs: differential expression genes. base.mbc.nctu.edu.tw/, accessed on 1 September 2020), miRDB [25] (http://www.mirdb.org/, accessed on 1 September 2020) and TargetScan [26] (http://www.targetscan.org/), the DEmRNA-miRNA interaction pairs were obtained. Only the overlapping of the interaction pairs supported by all databases was retained to create the miRNA-mRNA pairs. Subsequently, the correlation between DElncRNAs and DEmRNAs was analyzed by Pearson’s correlation coefficient (PCC) with the threshold of the PCC value r > 0.4 and p < 0.05 indicating a strong correlation. Finally, the remaining Biology 2021, 10, 432 5 of 15 A Limma package [21] (version 3.44.3; https://bioconductor.org/packages/release/ bioc/html/limma.html, accessed on 29 August 2020) in R was applied to identify the differentially expressed lncRNAs (DElncRNAs) and differentially expressed mRNAs (DEmRNAs) before a quantile normalization and base-2 logarithmic (log ) transforma- tion. The criteria were selected as a false discovery rate (FDR) < 0.05 and a fold change (FC) > 1.5. Volcano plots were performed for the DEGs using the ggplot2 package [22] (version 3.3.2; https://cran.r-project.org/web/packages/ggplot2/index.html, accessed on 29 August 2020) in the R software. Pheatmap package (version: 1.0.12; https://cran.r- project.org/web/packages/pheatmap/, accessed 29 August 2020) in R was used to gen- erate the hierarchical cluster heatmap representing the expression intensity and direction (Figure 1b). 2.3. Construction of the cSCC-Associated lncRNA-miRNA-mRNA Network Due to the absence of miRNA information, the miRcode database [23] (http://www. mircode.org/, accessed on 1 September 2020) was used to predict the targeted miR- NAs of DElncRNAs. According to the miRTarBase [24] (http://mirtarbase.mbc.nctu. edu.tw/, accessed on 1 September 2020), miRDB [25] (http://www.mirdb.org/, accessed on 1 September 2020) and TargetScan [26] (http://www.targetscan.org/), the DEmRNA- miRNA interaction pairs were obtained. Only the overlapping of the interaction pairs supported by all databases was retained to create the miRNA-mRNA pairs. Subsequently, the correlation between DElncRNAs and DEmRNAs was analyzed by Pearson’s correlation coefficient (PCC) with the threshold of the PCC value r > 0.4 and p < 0.05 indicating a strong correlation. Finally, the remaining miRNA-DElncRNA and miRNA-DEmRNA and their shared miRNAs were incorporated into the ceRNA network. Cytoscape [27] (version 3.6.1; http://cytoscape.org, accessed on 6 September 2020) was utilized to visualize the network. 2.4. Functional Enrichment Analysis The Gene Ontology (GO) functional enrichment analysis in the category Biologi- cal Processes (BP) and the KEGG pathway enrichment analysis of DEmRNAs in the ceRNA network were determined by the clusterProfiler package in R [28] (version 3.14.0; https://bioconductor.org/packages/release/bioc/html/clusterProfiler.html, accessed on 12 September 2020) and R. The threshold was selected to be p < 0.05. 2.5. Establishment of the Protein-Protein Interaction (PPI) Network and the Core ceRNA Subnetwork The Search Tool for the Retrieval of Interacting Genes (STRING) [29] (http://string- db.org/, accessed on 12 September 2020) was utilized to analyze PPI information with the threshold of a combined score >0.4. The plug-in Molecular Complex Detection (MCODE; version:1.4.2; http://apps.cytoscape.org/apps/mcode, accessed on 12 September 2020) [30] in Cytoscape was utilized to identify functionally-related hub modules with a degree cut- off of more than 5, a node score cut-off of 0.2, a k-score of 2 and maximum depth of 100. The key genes and related lncRNAs and miRNAs were then mapped into the core ceRNA subnetwork. 2.6. Expression Validation in an External Dataset To validate the expression of the DElncRNAs and DEmRNAs in the subnetwork, the GSE7553 based on the GPL570 platform with 11 SCC and 4 NS samples was obtained from GEO. The inclusion criteria of the cSCC samples, the data processing of the identifying DElncRNAs and DEmRNAs of GSE7553 datasets and the cut-off criteria were consistent with those of the GSE42677 and GSE45164 datasets. If the relative expression level of genes in the core ceRNA network was significantly different and met the same expression patterns with the GSE42677 and GSE45164 datasets, we considered those genes validated. The comparison of the log transformed expression level of the DElncRNAs and DEmRNAs was performed by GraphPad Prism V8.3.0 software (GraphPad Software, Inc., San Diego, CA, USA) with an unpaired t-test and p < 0.05 indicating a significant difference. Biology 2021, 10, 432 6 of 15 3. Results 3.1. A Total Of 24 DElncRNAs and 3221 DEmRNAs Were Identified between cSCC and NS Based on the PCA results of all samples from the GSE42677 and GSE45164 datasets, one cSCC sample and one NS sample affecting the classification were excluded for the down- stream analyses (Figure S1c). The heatmap of the DEGs showed that the gene expression levels between the remaining 19 cSCC and 12 NS samples were completely distinguishable (Figure 1c). Compared with the NS samples, 24 DElncRNAs (18 over-expressed and 6 down-expressed) and 3221 DEmRNAs (1568 over-expressed and 1653 down-expressed) were identified with the threshold of FDR < 0.05 and |FC| > 1.5. 3.2. The cSCC-Associated ceRNA Network Establishment Before predicting the miRNAs interacting with DElncRNAs or DEmRNAs, the co- expression analysis by the PCC was applied to calculate the correlation between DElncRNA and DEmRNA. Based on the result of the PCC with a cut-off value r > 0.4 and p < 0.05, 155 DEmRNAs and 23 DElncRNAs were retained for a further analysis (Table S1). The re- maining DEmRNAs and DElncRNAs that shared the common miRNAs were incorporated Biology 2021, 10, x into the construction of the ceRNA network (Figure 2). Finally, a cSCC-associated 7 ceRNA of 16 network containing 137 miRNA-lncRNA (Table S2) and 221 miRNA-mRNA (Table S3) pairs was created composed of 31 miRNAs, 11 DElncRNAs (8 down-expressed and 3 over- expressed, Table S4) and 155 DEmRNAs (130 down-expressed and 25 over-expressed) and 3 over-expressed, Table S4) and 155 DEmRNAs (130 down-expressed and 25 over- (Tables S4 and S5). expressed) (Tables S4 and S5). Figure 2. The construction of the ceRNA network including 11 DElncRNAs, 31 predicted miRNAs and 155 DEmRNAs. The Figure 2. The construction of the ceRNA network including 11 DElncRNAs, 31 predicted miRNAs and 155 DEmRNAs. The circle in green, the rectangle in red and the triangle in blue represent DEmRNAs, miRNAs and DElncRNAs, respectively. circle in green, the rectangle in red and the triangle in blue represent DEmRNAs, miRNAs and DElncRNAs, respectively. DElncRNAs: differentially expressed lncRNAs; DEmRNAs: differentially expressed mRNAs; miRNAs: microRNAs. DElncRNAs: differentially expressed lncRNAs; DEmRNAs: differentially expressed mRNAs; miRNAs: microRNAs. 3.3. DEmRNAs of the ceRNA Network Were Associated with Enhancer-Binding, Protein Kinase 3.3. DEmRNAs of the ceRNA Network Were Associated with Enhancer-Binding, Protein Kinase Regulator Activity and Autophagy Regulator Activity and Autophagy To uncover the biological mechanisms associated with the development of a cSCC, a To uncover the biological mechanisms associated with the development of a cSCC, GO functional and a KEGG pathway enrichment analysis of the DEmRNAs in the ceRNA a GO functional and a KEGG pathway enrichment analysis of the DEmRNAs in the network were performed. In the GO enrichment analysis (Figure 3a, Table S6), the DEmR- ceRNA network were performed. In the GO enrichment analysis (Figure 3a, Table S6), the NAs were involved in the biological processes of enhancer sequence-specific DNA bind- DEmRNAs were involved in the biological processes of enhancer sequence-specific DNA ing, enhancer-binding, protein kinase regulator activity and transcription corepressor ac- binding, enhancer-binding, protein kinase regulator activity and transcription corepressor tivity. The significant KEGG pathway (Figure 3b, Table S7) enrichment included the fork- activity. The significant KEGG pathway (Figure 3b, Table S7) enrichment included the head box O (FoxO) signaling pathway, autophagy and cellular senescence. forkhead box O (FoxO) signaling pathway, autophagy and cellular senescence. Biology 2021, 10, 432 7 of 15 Biology 2021, 10, x 8 of 16 (a) (b) Figure 3. Bar chart of the GO and KEGG analysis of DEmRNAs in the ceRNA network. (a) Bar chart of the significant GO Figure 3. Bar chart of the GO and KEGG analysis of DEmRNAs in the ceRNA network. (a) Bar chart of the significant GO terms with DEmRNAs identified by the clusterProfiler package in R. The color of the bar represents the −log10 transformed terms with DEmRNAs identified by the clusterProfiler package in R. The color of the bar represents the log transformed p values of the GO terms and the abscissa corresponds with the number of genes in the GO terms; (b) bar chart of signifi- p values of the GO terms and the abscissa corresponds with the number of genes in the GO terms; (b) bar chart of significant cant KEGG enriched pathways with DEmRNAs identified by the clusterProfiler package in R. The color of the bar repre- KEGG enriched pathways with DEmRNAs identified by the clusterProfiler package in R. The color of the bar represents sents the −log10 transformed p values of the KEGG enriched pathways and the abscissa corresponds with the number of the log10 transformed p values of the KEGG enriched pathways and the abscissa corresponds with the number of genes in the KEGG enriched pathways. GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; DElncRNAs: differentially expressed lncRNAs; DEmRNAs: differentially expressed mRNAs. Biology 2021, 10, x 9 of 16 Biology 2021, 10, 432 8 of 15 genes in the KEGG enriched pathways. GO: Gene Ontology; KEGG: Kyoto Encyclopedia of Genes and Genomes; DElncRNAs: differentially expressed lncRNAs; DEmRNAs: differentially expressed mRNAs. 3.4. PPI Network and the Core ceRNA Subnetwork Construction 3.4. PPI Network and the Core ceRNA Subnetwork Construction The combination of the STRING tool and Cytoscape software constructed the PPI The combination of the STRING tool and Cytoscape software constructed the PPI network. A total of 99 nodes and 170 interactions were screened from the PPI network. After network. A total of 99 nodes and 170 interactions were screened from the PPI network. a MCODE analysis, five hub-genes were incorporated into the PPI network (Figure 4a,b). After a MCODE analysis, five hub-genes were incorporated into the PPI network (Figure Based on the ceRNA network, the hub-genes associated with the core ceRNA subnetwork 4a,b). Based on the ceRNA network, the hub-genes associated with the core ceRNA sub- (Tables S8 and S9) were established including 20 miRNA-lncRNA and 7 miRNA-mRNA network (Tables S8 and S9) were established including 20 miRNA-lncRNA and 7 miRNA- edges and 19 nodes (9 lncRNAs (HCG9, POLR2J4, HLA-F-AS1, HCG18, DLEU2, AC074212, mRNA edges and 19 nodes (9 lncRNAs (HCG9, POLR2J4, HLA-F-AS1, HCG18, DLEU2, SNX29P2, LINC00115 and CTC-338M12), 5 miRNAs (hsa-mir-1297, hsa-mir-137, hsa-mir- AC074212, SNX29P2, LINC00115 and CTC-338M12), 5 miRNAs (hsa-mir-1297, hsa-mir- 17-5p, 137, hsa-mir hsa-mir -20b-5p -17-5p, hs and a-mir hsa-mir -20b-5p -25b-5p) and hsand a-mir 5-25b mRNAs -5p) an (AGO4, d 5 mRNA E2F1, s (AG CCND1, O4, E2 TNRC6B F1, CCND1, TNRC6B and ETS1)) (Figure 5). and ETS1)) (Figure 5). (a) (b) Figure 4. PPI network of DEmRNAs and Cluster 1 using the MCODE plug-in from the whole PPI network. (a) PPI network Figure 4. PPI network of DEmRNAs and Cluster 1 using the MCODE plug-in from the whole PPI network. (a) PPI network of DEmRNAs in the ceRNA network obtained from the STRING database. (b) The first cluster is of PPI followed by the of DEmRNAs in the ceRNA network obtained from the STRING database. (b) The first cluster is of PPI followed by the Biology 2021, 10, x 10 of 16 MCODE plug-in of Cytoscape. PPI: Protein-protein interaction; DEmRNAs: differentially expressed mRNAs; MCODE: MCODE plug-in of Cytoscape. PPI: Protein-protein interaction; DEmRNAs: differentially expressed mRNAs; MCODE: molecular complex. molecular complex. Figure 5. The hub-genes associated with the ceRNA network in a cSCC. MRNAs, miRNAs and Figure 5. The hub-genes associated with the ceRNA network in a cSCC. MRNAs, miRNAs and lncRNAs are indicated by the circle in green, rectangle in red and triangle in blue, respectively. lncRNAs are indicated by the circle in green, rectangle in red and triangle in blue, respectively. 3.5. Verification of the DEGs in the Core ceRNA Subnetwork The expression of 9 DElncRNAs and 5 DEmRNAs in the hub-genes associated with the core ceRNA subnetwork were compared between the cSCC and the NS group using the GSE7553 dataset. Finally, one lncRNA HLA-F antisense RNA 1 (HLA-F-AS1) and 3 DEmRNAs (Argonaute-4 (AGO4), Cyclin D1 (CCND1) and E2F transcription factor 1 (E2F1)) were considered to be validated of which the expression patterns were consistent with the results of GSE42677 and GSE45164 (Figure 6). (a) (b) Biology 2021, 10, x 10 of 16 Biology 2021, 10, 432 9 of 15 Figure 5. The hub-genes associated with the ceRNA network in a cSCC. MRNAs, miRNAs and lncRNAs are indicated by the circle in green, rectangle in red and triangle in blue, respectively. 3.5. Verification of the DEGs in the Core ceRNA Subnetwork 3.5. Verification of the DEGs in the Core ceRNA Subnetwork The expression of 9 DElncRNAs and 5 DEmRNAs in the hub-genes associated with The expression of 9 DElncRNAs and 5 DEmRNAs in the hub-genes associated with the core ceRNA subnetwork were compared between the cSCC and the NS group using the core ceRNA subnetwork were compared between the cSCC and the NS group using the GSE7553 dataset. Finally, one lncRNA HLA-F antisense RNA 1 (HLA-F-AS1) and 3 the GSE7553 dataset. Finally, one lncRNA HLA-F antisense RNA 1 (HLA-F-AS1) and 3 DEmRNAsDEm (Argonaute-4 RNAs (Argon (AGO4), auteCyclin -4 (AGO4 D1 (CCND1) ), Cyclin and D1 (CCN E2F transcription D1) and E2F factor transcription 1 factor 1 (E2F1)) were considered to be validated of which the expression patterns were consistent (E2F1)) were considered to be validated of which the expression patterns were consistent with the results of GSE42677 and GSE45164 (Figure 6). with the results of GSE42677 and GSE45164 (Figure 6). Biology 2021, 10, x 11 of 16 (a) (b) (c) (d) Figure 6. External validation of the four core genes in the GSE7553 datasets. (a) The relative Figure 6. External validation of the four core genes in the GSE7553 datasets. (a) The relative expression level of HLA-F- AS1 between expr cSCC ession and level NS sam of HLA-F-AS1 ples; (b) the between relativecSCC expres and sion NSle samples; vel of AGO4 (b) thebetwee relative n expr cSCC ession and level NS samples; of (c) the relative expression level of CCND1 between cSCC and NS samples; (d) the relative expression level of E2F1 between cSCC AGO4 between cSCC and NS samples; (c) the relative expression level of CCND1 between cSCC and NS samples. The log2 transformed expression level corresponded with the relative expression level. Data are presented and NS samples; (d) the relative expression level of E2F1 between cSCC and NS samples. The log as a median with interquartile ranges. An unpaired t-test was used for the comparison of the relative expression level transformed expression level corresponded with the relative expression level. Data are presented as between the two groups. a median with interquartile ranges. An unpaired t-test was used for the comparison of the relative expression level between the two groups. 4. Discussion A cutaneous squamous cell carcinoma (cSCC) derived from the malignant prolifera- tion of keratinocytes is the second most frequent cause of non-melanoma skin cancer [31]. However, a cSCC usually displays a primary cSCC with benign clinical behavior and can be successfully treated with surgery. It may progress to be locally invasive and metastatic [32]. The etiology and pathogenesis of a cSCC have not been well elucidated. Recently, accumulating studies have emphasized the regulatory role of lncRNAs as ceRNAs in the occurrence and progression of many types of cancers including cSCCs [15–18]. We then constructed an lncRNA-associated ceRNA network based on the transcriptome expres- sion profiles from the GEO database to identify new targets with a potential diagnostic or curative value for the cSCC. In total, 24 DElncRNAs and 3221 DEmRNAs were identified between the cSCC tu- mors compared with NS samples. The DElncRNAs and DEmRNAs with a strong relation- ship were incorporated with their common targeted miRNAs to construct the ceRNA net- work. The DEmRNAs in the ceRNA network were then subjected to the GO enrichment and KEGG pathway enrichment analyses. The pathway analysis showed that the AMPK signaling pathway [33] and autophagy [34] participate in the pathogenesis of a cSCC. We established a hub-gene-associated core ceRNA subnetwork including nine lncRNAs, five predicted miRNAs (hsa-mir-1297, hsa-mir-137, hsa-mir-17-5p, hsa-mir-20b- 5p and hsa-mir-25b-5p) and five mRNAs (AGO4, E2F1, CCND1, TNRC6B and ETS1). Three DEmRNAs (E2F1, CCND1 and AGO4) were validated in the independent dataset GSE7553. E2F1 and CCND1 are both cell cycle regulatory proteins. The upregulated ex- pression of E2F1 and its family number E2F7 in cSCCs may contribute to cSCC formation [35]. E2F1 was shown to inhibit the differentiation of keratinocytes at the terminal phase and inactivate MMP [36,37]. The aberrant activation of E2F1 also participated in the cellu- lar proliferation, differentiation and apoptosis of colon cancer and prostate cancer and is often related to a poor prognosis [38,39]. CCND1 has been discovered to be upregulated in SCCs such as laryngeal squamous cell carcinomas (LSCCs), head and neck squamous cell carcinomas (HNSCCs) and oral squamous cell carcinomas (OSCCs) [40–42]. As for the AGO4, it may exert a function as an effector protein with DNA methylation via RNA- Biology 2021, 10, 432 10 of 15 4. Discussion A cutaneous squamous cell carcinoma (cSCC) derived from the malignant prolifera- tion of keratinocytes is the second most frequent cause of non-melanoma skin cancer [31]. However, a cSCC usually displays a primary cSCC with benign clinical behavior and can be successfully treated with surgery. It may progress to be locally invasive and metastatic [32]. The etiology and pathogenesis of a cSCC have not been well elucidated. Recently, accumu- lating studies have emphasized the regulatory role of lncRNAs as ceRNAs in the occurrence and progression of many types of cancers including cSCCs [15–18]. We then constructed an lncRNA-associated ceRNA network based on the transcriptome expression profiles from the GEO database to identify new targets with a potential diagnostic or curative value for the cSCC. In total, 24 DElncRNAs and 3221 DEmRNAs were identified between the cSCC tumors compared with NS samples. The DElncRNAs and DEmRNAs with a strong relationship were incorporated with their common targeted miRNAs to construct the ceRNA network. The DEmRNAs in the ceRNA network were then subjected to the GO enrichment and KEGG pathway enrichment analyses. The pathway analysis showed that the AMPK signaling pathway [33] and autophagy [34] participate in the pathogenesis of a cSCC. We established a hub-gene-associated core ceRNA subnetwork including nine lncR- NAs, five predicted miRNAs (hsa-mir-1297, hsa-mir-137, hsa-mir-17-5p, hsa-mir-20b-5p and hsa-mir-25b-5p) and five mRNAs (AGO4, E2F1, CCND1, TNRC6B and ETS1). Three DEmRNAs (E2F1, CCND1 and AGO4) were validated in the independent dataset GSE7553. E2F1 and CCND1 are both cell cycle regulatory proteins. The upregulated expression of E2F1 and its family number E2F7 in cSCCs may contribute to cSCC formation [35]. E2F1 was shown to inhibit the differentiation of keratinocytes at the terminal phase and inactivate MMP [36,37]. The aberrant activation of E2F1 also participated in the cellular proliferation, differentiation and apoptosis of colon cancer and prostate cancer and is often related to a poor prognosis [38,39]. CCND1 has been discovered to be upregulated in SCCs such as laryngeal squamous cell carcinomas (LSCCs), head and neck squamous cell carcinomas (HNSCCs) and oral squamous cell carcinomas (OSCCs) [40–42]. As for the AGO4, it may exert a function as an effector protein with DNA methylation via RNA-dependent DNA methylation (RdDM) [43]. A DNA methyltransferase 3A (DNMT3A)/AGO4 complex abol- ished miR-181-5p inhibition of gene expression via the cytosine methylation of miR-181-5p and resulted in an aggressive outcome of a glioblastoma multiforme (GBM) [44]. Among the five predicted miRNAs, miR-17-5p had been validated to be highly and significantly upregulated in a cSCC by a real-time polymerase chain reaction (RT-PCR) between cSCC tissues and NS samples [45]. MiR-17-5p and miR-20b-5p are both members of the miR-17-92 family. Moreover, the interaction of miR-17 with CCND1 and E2F1 has been studied. In a study of a lung adenocarcinoma, miR-20b-5p was found to directly and independently inhibit the activity of reporter genes and decrease the CCND1 protein 0 0 expression [46]. MiR-20b-5p could sponge with the 3 untranslated region (3 UTR) of CCND1, repressing the CCND1 expression in a lung cancer cell [46]. MiR-20b-5p also targeted with E2F1, inducing a myoblast differentiation and repressing its proliferation [47]. Only one lncRNA HLA-F-AS1 in the core ceRNA subnetwork was validated in the external datasets. Several studies had revealed the dysregulation of HLA-F-AS1 in can- cers [48–50] while the transcript expression of HLA-F-AS1 was verified to be downregulated in lung adenocarcinoma (LAD) samples [49]. It was over-expressed in colorectal cancer (CRC) tissues and CRC cell lines. Huang et al. found that HLA-F-AS1 may participate in the progression of CRC via the miR-330-3p/profilin 1 (PFN1) axis in CRC cells [48]. The elevated HLA-F-AS1 expression induced by STAT3 promoted cell growth and stemness characteristics via sponging with miR-541-3p to upregulate the expression of a TraB domain (TRABD) containing triple-negative breast cancer cells [50]. In summary, according to the predicted core ceRNA subnetwork, lncRNA HLA-F-AS1 may sponge miR-17-5p or miR-20b-5p to regulate the expression of CCND1 and E2F1 in the cSCC. The expression of HLA-F-AS1 in the cSCC samples and its role mechanism in Biology 2021, 10, 432 11 of 15 the biological process of a cSCC as a ceRNA required further research. To highlight the main findings of the validated lncRNAs and miRNAs in the core ceRNA subnetwork, we followed the model presented by Garofoli et al. [51] and made a table where a systematic review of the roles and mechanisms of HLA-F-AS1, miR-17-5p and miR-20b-5p in the pathogenesis of cancers was carried out (Table 1). Table 1. Summary of the roles and mechanisms of discussed lncRNAs and miRNAs in the pathogenesis of cancers. LncRNA/miRNA Gene Targets Mechanisms in Cancers Promotion of cell proliferation and stemness in vitro and tumor growth in vivo by MiR-541-3p the miR-541-3p/TRABD axis in triple-negative breast cancer [50]. MiR-330-3p HLA-F-AS1 Promotion of cell proliferation, migration and invasion and repression of apoptosis in vitro via the miR-330-3p/ PFN1 axis in colorectal cancer [48]. 0 0 3 UTR of CCND1 Repression of the CCND1 protein expression via sponging with the CCND1 3 -UTR E2F1 region in a lung cancer cell [46]. MiR-17-5p Promotion of cell growth, proliferation, migration and colony formation in vitro in a non-metastatic hepatocellular carcinoma by regulating E2F1 expression [52]. Promotion of myoblast differentiation and repression of myoblast proliferation by 3 UTR of CCND1 regulating E2F1 expression [47]. 3 UTR of E2F1 Promotion of cell proliferation, G1/S transition and DNA synthesis by downregulating CCND1 and E2F1 expression in multipotent stromal cells [53]. MiR-20b-5p Promotion of cell migration and proliferation and repression of apoptosis of breast cancer stem cells and promoting tumor growth in vivo by bidirectionally regulating the protein levels of CCND1 and E2F1 in breast cancer [54]. Inhibition of the cell cycle, migration and invasion in vitro and the tumorigenesis in vivo by negatively regulating CCND1 in colon cancer [55]. HLA-F-AS1: HLA-F antisense RNA 1; CCND1: Cyclin D1; E2F1: E2F transcription factor 1; TRABD: TraB domain containing; PFN1: 0 0 profilin 1; 3 UTR: 3 untranslated regions. 5. Conclusions In this study, the cSCC-related lncRNA-miRNA-mRNA ceRNA network was es- tablished and one lncRNA (HLA-F-AS1) and three mRNAs (AGO4, CCND1 and E2F1) in the core ceRNA network were validated. The HLA-F-AS1-miR-17-5p/miR-20b-5p- CCND1/E2F1 axis may take part in the development of a cSCC. However, there are several limitations in our research that should be considered. First, miRNA expression datasets are needed to verify the network. Second, a relatively large sample size should be collected to assure the results are more reliable in future studies. Moreover, the relation between the ceRNA network and clinical features could not be established due to the absence of related clinical information. Our research provides a new therapeutic idea and regulatory roles of ceRNAs in the pathogenesis of a cSCC. They may be candidate targets with a promising diagnostic or therapeutic value with a cSCC. Further functional investigations mediated by the ceRNA would have to be verified. Supplementary Materials: The following are available online at https://www.mdpi.com/article/ 10.3390/biology10050432/s1, Figure S1: Data preprocessing. Table S1: The lncRNA-mRNA pairs with significant correlations, Table S2: DElncRNA-miRNA pairs of the ceRNA network, Table S3: DEmRNA-miRNA pairs of the ceRNA network, Table S4: Over-expressed and down-expressed lncR- NAs of the ceRNA network, Table S5: Over-expressed and down-expressed mRNAs of the ceRNA network, Table S6: GO enrichment terms of DEmRNAs in the ceRNA network, Table S7: KEGG pathway enrichment of DEmRNAs in the ceRNA network, Table S8: Details of DElncRNA-miRNA pairs of the core ceRNA network, Table S9: DEmRNA-miRNA pairs of the core ceRNA network. Author Contributions: Y.X., Y.D. (Yunhua Deng) and Q.Q. conceived the idea. Y.X., Y.D. (Yingying Dong), M.W., H.L., Q.S., Q.G. and Q.Q. performed the literature review, data analysis and wrote the manuscript. Y.D. (Yunhua Deng) and Q.Q. revised the manuscript and helped to improve the quality of the manuscript. All authors have read and agreed to the published version of the manuscript. Biology 2021, 10, 432 12 of 15 Funding: This research was funded by the National Natural Science Foundation of China (Grant No. 81773306). Institutional Review Board Statement: Not applicable. Informed Consent Statement: Not applicable. Data Availability Statement: Data were obtained from the GEO database and are available at: www.ncbi.nlm.nih.gov/geo (accessed on 24 August 2020) with the permission of the GEO database. Conflicts of Interest: The authors declare no conflict of interest. Abbreviations SCC Squamous Cell Carcinoma cSCC Cutaneous Squamous Cell Carcinoma SCCIS Squamous Cell Carcinoma (SCC) In Situ NcRNAs Non-Coding RNAs LncRNAs Long Non-Coding RNAs MiRNAs MicroRNAs CeRNAs Competing Endogenous RNAs NS Normal Skin DEGs Differentially Expressed Genes DElncRNAs Differentially Expressed LncRNAs DEmRNAs Differentially Expressed MRNAs GEO Gene Expression Omnibus GO Gene Ontology BP Biological Processes KEGG Kyoto Encyclopedia of Genes and Genomes STRING Search Tool for the Retrieval of Interacting Genes MCODE Molecular Complex Detection FDR False Discovery Rate FC Fold Change PCA Principal Component Analysis PCC Pearson’s Correlation Coefficient CDC20 Cell Division Cycle 20 EGFR Epidermal Growth Factor Receptor PRAF2 Prenylated Rab Acceptor 1 Domain Family Member 2 FoxO Forkhead box O HLA-F-AS1 HLA-F Antisense RNA 1 CCND1 Cyclin D1 E2F1 E2F Transcription Factor 1 AGO4 Argonaute-4 LAD Lung Adenocarcinoma CRC Colorectal Cancer PFN1 Profilin 1 ESCC Esophageal Squamous Cell Carcinoma TSCC Tongue Squamous Cell Carcinoma Progression LSCC Laryngeal Squamous Cell Carcinoma HNSCC Head and Neck Squamous Cell Carcinoma OSCC Oral Squamous Cell Carcinoma RdDM RNA-Dependent DNA Methylation DNMT3A DNA Methyltransferase 3A GBM Glioblastoma Multiforme AK Actinic Keratosis NEAT1 Nuclear Paraspeckle Assembly Transcript 1 RT-PCR Real-time Polymerase Chain Reaction 0 0 3 UTRs 3 Untranslated Regions TRABD TraB Domain Containing Biology 2021, 10, 432 13 of 15 References 1. Motaparthi, K.; Kapil, J.P.; Velazquez, E.F. Cutaneous Squamous Cell Carcinoma: Review of the Eighth Edition of the American Joint Committee on Cancer Staging Guidelines, Prognostic Factors, and Histopathologic Variants. Adv. Anat. Pathol. 2017, 24, 171–194. [CrossRef] 2. Waldman, A.; Schmults, C. Cutaneous Squamous Cell Carcinoma. Hematol. Oncol. Clin. N. Am. 2019, 33, 1–12. [CrossRef] 3. Fernandez Figueras, M.T. From actinic keratosis to squamous cell carcinoma: Pathophysiology revisited. J. Eur. Acad. Derm. Venereol. 2017, 31, 5–7. [CrossRef] [PubMed] 4. Soura, E.; Gagari, E.; Stratigos, A. Advanced cutaneous squamous cell carcinoma: How is it defined and what new therapeutic approaches are available. Curr. Opin. Oncol. 2019, 31, 461–468. [CrossRef] [PubMed] 5. Esteller, M. Non-coding RNAs in human disease. Nat. Rev. Genet. 2011, 12, 861–874. [CrossRef] 6. Bai, X.; Zhou, Y.; Chen, P.; Yang, M.; Xu, J. MicroRNA-142-5p induces cancer stem cell-like properties of cutaneous squamous cell carcinoma via inhibiting PTEN. J. Cell. Biochem. 2018, 119, 2179–2188. [CrossRef] [PubMed] 7. Wang, Y.; Deng, X.; Dai, Y.; Niu, X.; Zhou, M. miR-27a Downregulation Promotes Cutaneous Squamous Cell Carcinoma Progression via Targeting EGFR. Front. Oncol. 2020, 9, 1565. [CrossRef] 8. Quinn, J.J.; Chang, H.Y. Unique features of long non-coding RNA biogenesis and function. Nat. Rev. Genet. 2016, 17, 47–62. [CrossRef] 9. Rotondo, J.C.; Selvatici, R.; Di Domenico, M.; Marci, R.; Vesce, F.; Tognon, M.; Martini, F. Methylation loss at H19 imprinted gene correlates with methylenetetrahydrofolate reductase gene promoter hypermethylation in semen samples from infertile males. Epigenetics 2013, 8, 990–997. [CrossRef] 10. Di Martino, M.T.; Riillo, C.; Scionti, F.; Grillone, K.; Polerà, N.; Caracciolo, D.; Arbitrio, M.; Tagliaferri, P.; Tassone, P. miRNAs and lncRNAs as Novel Therapeutic Targets to Improve Cancer Immunotherapy. Cancers (Basel) 2021, 13, 1587. [CrossRef] 11. Mei, X.L.; Zhong, S. Long noncoding RNA LINC00520 prevents the progression of cutaneous squamous cell carcinoma through the inactivation of the PI3K/Akt signaling pathway by downregulating EGFR. Chin. Med. J. (Engl.) 2019, 132, 454–465. [CrossRef] 12. Piipponen, M.; Heino, J.; Kähäri, V.M.; Nissinen, L. Long non-coding RNA PICSAR decreases adhesion and promotes migration of squamous carcinoma cells by downregulating 2 1 and 5 1 integrin expression. Biol. Open 2018, 7, bio037044. [CrossRef] 13. Yamamura, S.; Imai-Sumida, M.; Tanaka, Y.; Dahiya, R. Interaction and cross-talk between non-coding RNAs. Cell. Mol. Life Sci. 2018, 75, 467–484. [CrossRef] 14. López-Urrutia, E.; Bustamante Montes, L.P.; Ladrón de Guevara Cervantes, D.; Pérez-Plasencia, C.; Campos-Parra, A.D. Crosstalk Between Long Non-coding RNAs, Micro-RNAs and mRNAs: Deciphering Molecular Mechanisms of Master Regulators in Cancer. Front. Oncol. 2019, 9, 669. [CrossRef] 15. Li, Y.; Chen, D.; Gao, X.; Li, X.; Shi, G. LncRNA NEAT1 Regulates Cell Viability and Invasion in Esophageal Squamous Cell Carcinoma through the miR-129/CTBP2 Axis. Dis. Markers 2017, 9, 5314649. [CrossRef] 16. Ding, J.; Yang, C.; Yang, S. LINC00511 interacts with miR-765 and modulates tongue squamous cell carcinoma progression by targeting LAMC2. J. Oral. Pathol. Med. 2018, 47, 468–476. [CrossRef] [PubMed] 17. Zhao, J.; Bai, X.; Feng, C.; Shang, X.; Xi, Y. Long Non-Coding RNA HCP5 Facilitates Cell Invasion And Epithelial-Mesenchymal Transition In Oral Squamous Cell Carcinoma By miR-140-5p/SOX4 Axis. Cancer Manag. Res. 2019, 11, 10455–10462. [CrossRef] [PubMed] 18. Yu, G.J.; Sun, Y.; Zhang, D.W.; Zhang, P. Long non-coding RNA HOTAIR functions as a competitive endogenous RNA to regulate PRAF2 expression by sponging miR-326 in cutaneous squamous cell carcinoma. Cancer Cell Int. 2019, 19, 270. [CrossRef] 19. Johnson, W.E.; Li, C.; Rabinovic, A. Adjusting batch effects in microarray expression data using empirical Bayes methods. Biostatistics 2007, 8, 118–127. [CrossRef] [PubMed] 20. Konishi, T. Principal component analysis for designed experiments. BMC Bioinform. 2015, 16, S7. [CrossRef] 21. Ritchie, M.E.; Phipson, B.; Wu, D.; Hu, Y.; Law, C.W.; Shi, W.; Smyth, G.K. limma powers differential expression analyses for RNA-sequencing and microarray studies. Nucleic Acids Res. 2015, 43, e47. [CrossRef] [PubMed] 22. Wickham, H. ggplot2: Elegant Graphics for Data Analysis; Springer: New York, NY, USA, 2016. 23. Jeggari, A.; Marks, D.S.; Larsson, E. miRcode: A map of putative microRNA target sites in the long non-coding transcriptome. Bioinformatics 2012, 28, 2062–2063. [CrossRef] 24. Chou, C.H.; Shrestha, S.; Yang, C.D.; Chang, N.W.; Lin, Y.L.; Liao, K.W.; Huang, W.C.; Sun, T.H.; Tu, S.J.; Lee, W.H.; et al. miRTarBase update 2018: A resource for experimentally validated microRNA-target interactions. Nucleic Acids Res. 2018, 46, D296–D302. [CrossRef] [PubMed] 25. Chen, Y.; Wang, X. miRDB: An online database for prediction of functional microRNA targets. Nucleic Acids Res. 2020, 48, D127–D131. [CrossRef] [PubMed] 26. Lewis, B.P.; Shih, I.-H.; Jones-Rhoades, M.W.; Bartel, D.P.; Burge, C.B. Prediction of mammalian microRNA targets. Cell 2003, 115, 787–798. [CrossRef] 27. Shannon, P.; Markiel, A.; Ozier, O.; Baliga, N.S.; Wang, J.T.; Ramage, D.; Amin, N.; Schwikowski, B.; Ideker, T. Cytoscape: A software environment for integrated models of biomolecular interaction networks. Genome Res. 2003, 13, 2498–2504. [CrossRef] [PubMed] 28. Yu, G.; Wang, L.-G.; Han, Y.; He, Q.-Y. clusterProfiler: An R package for comparing biological themes among gene clusters. Omics 2012, 16, 284–287. [CrossRef] Biology 2021, 10, 432 14 of 15 29. Szklarczyk, D.; Franceschini, A.; Wyder, S.; Forslund, K.; Heller, D.; Huerta-Cepas, J.; Simonovic, M.; Roth, A.; Santos, A.; Tsafou, K.P.; et al. STRING v10: Protein-protein interaction networks, integrated over the tree of life. Nucleic Acids Res. 2015, 43, D447–D452. [CrossRef] 30. Bader, G.D.; Hogue, C.W.V. An automated method for finding molecular complexes in large protein interaction networks. BMC Bioinform. 2003, 4, 2. [CrossRef] 31. Que, S.K.T.; Zwald, F.O.; Schmults, C.D. Cutaneous squamous cell carcinoma: Management of advanced and high-stage tumors. J. Am. Acad. Derm. 2018, 78, 249–261. [CrossRef] 32. Brougham, N.D.; Dennett, E.R.; Cameron, R.; Tan, S.T. The incidence of metastasis from cutaneous squamous cell carcinoma and the impact of its risk factors. J. Surg. Oncol. 2012, 106, 811–815. [CrossRef] 33. Piipponen, M.; Nissinen, L.; Farshchian, M.; Riihilä, P.; Kivisaari, A.; Kallajoki, M.; Peltonen, J.; Peltonen, S.; Kähäri, V.M. Long Noncoding RNA PICSAR Promotes Growth of Cutaneous Squamous Cell Carcinoma by Regulating ERK1/2 Activity. J. Investig. Derm. 2016, 136, 1701–1710. [CrossRef] [PubMed] 34. Zhou, W.; Zhang, S.; Li, J.; Li, Z.; Wang, Y.; Li, X. lncRNA TINCR participates in ALA-PDT-induced apoptosis and autophagy in cutaneous squamous cell carcinoma. J. Cell. Biochem. 2019, 120, 13893–13902. [CrossRef] 35. Endo-Munoz, L.; Dahler, A.; Teakle, N.; Rickwood, D.; Hazar-Rethinam, M.; Abdul-Jabbar, I.; Sommerville, S.; Dickinson, I.; Kaur, P.; Paquet-Fifield, S.; et al. E2F7 can regulate proliferation, differentiation, and apoptotic responses in human keratinocytes: Implications for cutaneous squamous cell carcinoma formation. Cancer Res. 2009, 69, 1800–1808. [CrossRef] [PubMed] 36. Das Mahapatra, K.; Pasquali, L.; Søndergaard, J.N.; Lapins, J.; Nemeth, I.B.; Baltás, E.; Kemény, L.; Homey, B.; Moldovan, L.-I.; Kjems, J.; et al. A comprehensive analysis of coding and non-coding transcriptomic changes in cutaneous squamous cell carcinoma. Sci. Rep. 2020, 10, 3637. [CrossRef] [PubMed] 37. Nagarajan, P.; Chin, S.S.; Wang, D.; Liu, S.; Sinha, S.; Garrett-Sinha, L.A. Ets1 blocks terminal differentiation of keratinocytes and induces expression of matrix metalloproteases and innate immune mediators. J. Cell Sci. 2010, 123, 3566–3575. [CrossRef] 38. Chun, J.N.; Cho, M.; Park, S.; So, I.; Jeon, J.-H. The conflicting role of E2F1 in prostate cancer: A matter of cell context or interpretational flexibility? Biochim. Biophys. Acta Rev. Cancer 2020, 1873, 188336. [CrossRef] 39. Fang, Z.; Lin, M.; Li, C.; Liu, H.; Gong, C. A comprehensive review of the roles of E2F1 in colon cancer. Am. J. Cancer Res. 2020, 10, 757–768. 40. Hanken, H.; Gröbe, A.; Cachovan, G.; Smeets, R.; Simon, R.; Sauter, G.; Heiland, M.; Blessmann, M. CCND1 amplification and cyclin D1 immunohistochemical expression in head and neck squamous cell carcinomas. Clin. Oral Investig. 2014, 18, 269–276. [CrossRef] 41. Wang, X.; Wang, C.; Yan, G.; Kang, Y.; Sun, G.; Wang, S.; Zou, R.; Sun, H.; Zeng, K.; Song, H.; et al. BAP18 is involved in upregulation of CCND1/2 transcription to promote cell growth in oral squamous cell carcinoma. EBioMedicine 2020, 53, 102685. [CrossRef] 42. Zang, Y.; Li, J.; Wan, B.; Tai, Y. circRNA circ-CCND1 promotes the proliferation of laryngeal squamous cell carcinoma through elevating CCND1 expression via interacting with HuR and miR-646. J. Cell. Mol. Med. 2020, 24, 2423–2433. [CrossRef] 43. Chalertpet, K.; Pin-On, P.; Aporntewan, C.; Patchsung, M.; Ingrungruanglert, P.; Israsena, N.; Mutirangura, A. Argonaute 4 as an Effector Protein in RNA-Directed DNA Methylation in Human Cells. Front. Genet. 2019, 10, 645. [CrossRef] 44. Cheray, M.; Etcheverry, A.; Jacques, C.; Pacaud, R.; Bougras-Cartron, G.; Aubry, M.; Denoual, F.; Peterlongo, P.; Nadaradjane, A.; Briand, J.; et al. Cytosine methylation of mature microRNAs inhibits their functions and is associated with poor prognosis in glioblastoma multiforme. Mol. Cancer 2020, 19, 36. [CrossRef] 45. Sand, M.; Hessam, S.; Amur, S.; Skrygan, M.; Bromba, M.; Stockfleth, E.; Gambichler, T.; Bechara, F.G. Expression of oncogenic miR-17-92 and tumor suppressive miR-143-145 clusters in basal cell carcinoma and cutaneous squamous cell carcinoma. J. Derm. Sci. 2017, 86, 142–148. [CrossRef] 46. Ding, H.; Luo, Y.; Hu, K.; Liu, P.; Xiong, M. Linc00467 promotes lung adenocarcinoma proliferation via sponging miR-20b-5p to activate CCND1 expression. Onco Targets Ther. 2019, 12, 6733–6743. [CrossRef] 47. Luo, W.; Li, G.; Yi, Z.; Nie, Q.; Zhang, X. E2F1-miR-20a-5p/20b-5p auto-regulatory feedback loop involved in myoblast proliferation and differentiation. Sci. Rep. 2016, 6, 27904. [CrossRef] [PubMed] 48. Huang, Y.; Sun, H.; Ma, X.; Zeng, Y.; Pan, Y.; Yu, D.; Liu, Z.; Xiang, Y. HLA-F-AS1/miR-330-3p/PFN1 axis promotes colorectal cancer progression. Life Sci. 2020, 254, 117180. [CrossRef] 49. Yang, Z.; Li, H.; Wang, Z.; Yang, Y.; Niu, J.; Liu, Y.; Sun, Z.; Yin, C. Microarray expression profile of long non-coding RNAs in human lung adenocarcinoma. Thorac. Cancer 2018, 9, 1312–1322. [CrossRef] [PubMed] 50. Wu, D.; Jia, H.; Zhang, Z.; Li, S. STAT3-induced HLA-F-AS1 promotes cell proliferation and stemness characteristics in triple negative breast cancer cells by upregulating TRABD. Bioorg. Chem. 2021, 109, 104722. [CrossRef] [PubMed] 51. Garofoli, M.; Volpicella, M.; Guida, M.; Porcelli, L.; Azzariti, A. The Role of Non-Coding RNAs as Prognostic Factor, Predictor of Drug Response or Resistance and Pharmacological Targets, in the Cutaneous Squamous Cell Carcinoma. Cancers (Basel) 2020, 12, 2552. [CrossRef] 52. El Tayebi, H.M.; Omar, K.; Hegy, S.; El Maghrabi, M.; El Brolosy, M.; Hosny, K.A.; Esmat, G.; Abdelaziz, A.I. Repression of miR-17-5p with elevated expression of E2F-1 and c-MYC in non-metastatic hepatocellular carcinoma and enhancement of cell growth upon reversing this expression pattern. Biochem. Biophys. Res. Commun. 2013, 434, 421–427. [CrossRef] [PubMed] Biology 2021, 10, 432 15 of 15 53. Tai, L.; Huang, C.J.; Choo, K.B.; Cheong, S.K.; Kamarul, T. Oxidative Stress Down-Regulates MiR-20b-5p, MiR-106a-5p and E2F1 Expression to Suppress the G1/S Transition of the Cell Cycle in Multipotent Stromal Cells. Int. J. Med. Sci. 2020, 17, 457–470. [CrossRef] [PubMed] 54. Xia, L.; Li, F.; Qiu, J.; Feng, Z.; Xu, Z.; Chen, Z.; Sun, J. Oncogenic miR-20b-5p contributes to malignant behaviors of breast cancer stem cells by bidirectionally regulating CCND1 and E2F1. BMC Cancer 2020, 20, 949. [CrossRef] 55. Yang, H.; Lin, J.; Jiang, J.; Ji, J.; Wang, C.; Zhang, J. miR-20b-5p functions as tumor suppressor microRNA by targeting cyclinD1 in colon cancer. Cell Cycle 2020, 19, 2939–2954. [CrossRef] [PubMed]

Journal

BiologyMultidisciplinary Digital Publishing Institute

Published: May 13, 2021

Keywords: competing endogenous RNAs network; long non-coding RNAs; microRNAs; cutaneous squamous cell carcinoma

There are no references for this article.