Oral squamous cell carcinoma (OSCC) is a well-known malignancy that frequently arises from the epidermal layer in oral cavity and, in most of the cases, develops from precancerous lesions of the oral mucosa. The majority of the OSCC patients are among those affected by head and neck squamous cell carcinoma (HNSCC), with an overall 5-year survival rate of ~64.4%. The outcome of OSCC treatment may vary depending on a series of factors, such as age, race, tumor staging, secondary complications and location of the tumor in the oral cavity (Zanoni et al., 2019). So far, current therapeutic strategies have been unable to fundamentally predict and fully cure OSCC patients. The mechanisms involved in the transformation of the normal oral epithelium into OSCC have not been fully characterized (Abdalla et al., 2017), therefore limiting any prediction on tumor progression and potential therapeutic approaches. Despite the various methods used for identifying the tumor type/stage and corresponding prognosis, most of these approaches have limited use to precisely specify the occurrence and development of the disease. Due to postoperative recurrence and lack of targeted chemical drugs, further innovations are urgently needed in regard to OSCC diagnosis and treatment. Therefore, a better prediction method for OSCC prognosis has been proposed, namely using gene expression data combined with clinical information of affected patients. Since microarray analyses and high-throughput sequencing techniques have been largely utilized in oncology, a representative number of genetic and molecular alterations have been reported (Shridhar et al., 2016). Still, the use of single markers to accurately diagnose malignant tumors is challenging. Thus, it appears critical to collectively distinguish candidate genes that are essential in cancer progression, from unrelated ones.
In general terms, Bioinformatics refers to a cutting-edge method of biological information analysis with obvious advantages when compared with traditional trials, such as the high efficiency to analyze extensive amounts of supporting data. As a result, Bioinformatic approaches have been continuously applied in cancer research to eliminate inconsistent results led by small sample size and/or by applying different technological platforms. However, the application of Bioinformatic techniques in the field of OSCC research is still restricted, despite the multitude of available information awaiting to be mined in public databases. Instead of performing isolated bioinformatics analysis, we have presently added biopsy RNA sequencing from OSCC patients during the analytical process. We further confronted this analysis with large sampling of microarray data in the current databases, accompanied with immunohistochemical (IHC) staining, therefore increasing the reliability of endpoint results. According to this Bioinformatic workflow, Small Proline Rich Protein 3 (SPRR3) was screened out as a potential prognostic and diagnostic biomarker for OSCC.
Small Proline Rich Protein 3 is a protein coding gene largely expressed in normal epithelial cells, whose gene product is known as a small proline-rich protein typically present in epidermal differentiation complexes. Importantly, SPRR3 has been recognized as a marker of terminal squamous cell differentiation (Candi, Schmidt & Melino, 2005), but its role in OSCC has been rarely studied. In our current work, we have established a bioinformatic-based workflow (Fig. 1) to screen novel diagnostic and prognostic OSCC markers as well as to verify their predictive efficiency in OSCC patients. As a result, we were able to identify SPRR3 as a putative target marker in OSCC.
Materials and Method
Information of specimens
We collected three pairs of tumor and adjacent oral normal specimens from three patients diagnosed with OSCC, whose diagnoses were supported by pathological examination results. All the three patients, coming from Stomatological Hospital of Shandong University, were considered as candidates for curative surgical resection (patients’ details in Table S1). This study was approved by the Research Ethics Committee of Stomatological Hospital of Shandong University (No. 20190205). Tissue microarray chips consisted of 61 OSCC samples and 10 samples of normal control was constructed by Shanghai Qutdo Biotech Company (Shanghai, China. Patients’ details in Table S2). The tissue microarray construction was approved by Taizhou Hospital Research Ethics Committee in Zhejiang Province (YB M-05-01). All experiments were performed complying with the relevant regulations, and written informed consent was obtained from patients.
High-throughput RNA sequence
We performed poly-A sequencing on the three pairs of tissues. According to the manufacturer’s protocol, we extracted total RNA from the three pairs of tissues we collected mentioned above by Trizol (Invitrogen, Carlsbad, CA, USA), and removed ribosomal RNA by the Ribo-Zero™ kit (Epicentre, Madison, WI, USA). The first and second strand cDNA synthesis of the fragment RNA (the average length of 200 bp) was performed, followed by adaptor ligation, and low cycle enrichment in accordance with the instructions of NEBNext® Ultra™ RNA Library Prep Kit for Illumina (New England Biolabs, Ipswich, MA, USA). In addition, the clean data quality control was filtered using Trimmomatic (Version 3). We utilized the Agilent 2200 TapeStation (Life Technologies, Waltham, MA, USA) and Qubit®2.0 (Invitrogen, Carlsbad, CA, USA) to evaluate the purified library products. The libraries were paired-end sequenced (PE150, Sequencing reads were 150 bp) at Guangzhou RiboBio Co., Ltd. (Guangzhou, Guangdong province, China) using IlluminaHiSeq 3000 platform (NEB).
Paired-end reads alignment to the human reference genome hg19 (NCBI build 37.1, http://www.ncbi.nlm.nih.gov/projects/genome/assembly/grc/) were performed with HISAT2 (http://ccb.jhu.edu/software/hisat2/manual.shtml#getting-started-with-hisat2, Pertea et al., 2016). We use HTSeq v0.6.0 (Anders, Pyl & Huber, 2015) to count the reads numbers mapped to each gene. The expression levels were presented as RPKM (expected number of Reads PerKilobase of transcript sequence per Million base pairs sequenced), which was the recommended and most frequently-used measuring unit. And the RNA-seq data was uploaded to GEO (Gene Expression Omnibus) database (GSE140707).
Data processing for public databsets
Gene Expression Profiling Interactive Analysis (GEPIA) (Tang et al., 2017) is an interactive web server for analyzing the RNA-seq expression data from the TCGA (the Cancer Genome Atlas) (https://tcga-data.nci.nih.gov/tcga/) (Tomczak, Czerwińska & Wiznerowicz, 2015) and Genotype-Tissue Expression (GTEx) projects. GEPIA provides customizable functions such as tumor/normal differential expression analysis, profiling according to cancer types or pathological stages, patient survival analysis, etc.
GSE3524 (Toruner et al., 2004), GSE42743 (Lohavanichbutr et al., 2013) and GSE30784 (Chen et al., 2008) datasets in the GEO database were selected as data source for OSCC information. Among them, GSE30784 and GSE3524 are the datasets containing gene expression information in both tumor and normal oral tissues, and GSE42743 is the one containing patients’ survival and clinical data such as survival time, live status, and other related factors that may influence the development and occurrence of OSCC. All of the datasets were selected to increase the sample size of the RNA-seq microarray data, so as to better verify the differential expression and diagnose the prognosis effect of the candidate genes. The platform files and original files (.CEL files) were downloaded. With the aid of Robust Muliti-array Average (RMA) algorithm, data imputation, background correction, log2 conversion, quantile normalization, pro summarization and missing value supplementation were performed on the matrix data of each GEO dataset through the “impute” package and “affy” package (Gautier et al., 2004) in R/Bioconductor software (Szklarczyk et al., 2011) (version 3.5.3). A collection of 306 specimens of OSCC was obtained from TCGA database, the level 3 RNA sequence (RNA-seq) data (.count files) was downloaded from TCGA website (https://www.cancer.gov/about-nci/organization/ccg/research/structural-genomics/tcga), and their correlated clinical information from FireBrowse database (http://firebrowse.org/) (Chabanais et al., 2018). The RNA-seq data in the TCGA database was normalized and processed using the “edgR” package (Zhang et al., 2019) prior to log2 conversion (Yang et al., 2019).
Differential expression analysis
All the differentially expressed genes (DEGs) were obtained through high-throughput RNA-seq. The statistically significant DEGs were acquired by utilizing the “DEseq” package (Wang et al., 2010), and set adjusted P-value (FDR) < 0.05 and |log2 (fold change (FC)) | > 1 as inclusion criteria. Then a hierarchical clustering analysis was performed using the “gplots” package (Ma et al., 2017) in R/Bioconductor software based on the expression level of DEGs in different groups, and colors represent different clustering information.
Construction of a protein–protein interaction network and module analysis
To assess the inter-relationships among DEGs, at protein level, these were mapped using the Search Tool for the Retrieval of Interacting Genes (STRING) database (Szklarczyk et al., 2017). The inter-relation among DEGs (combined score ≥ 0.4) was utilized to design a protein–protein interaction (PPI) network using Cytoscape (Shannon et al., 2003) (v3.6.1). The molecular complex detection (MCODE) plugin was used to select meaningful modules in the PPI network. For this, MCODE k-core > 2 and MCODE score > 5 were utilized as cut-off values. MCODE is a widely used algorithm for predicting molecular complexes, according to the inter-connectivity and density between nodes in PPI networks. The respective score determines the number of included nodes in clusters as well as the cluster size. An appropriate score is more likely to help select the meaningful cluster, which bears potential biological functions. The k-core is a vertex weighting method. The weighting scheme further boosts the weight of densely connected vertices, to further optimize the module-based cluster scoring system. A higher k-core of a particular network is the central and most densely connected subnet (Bader & Hogue, 2003).
Kyoto encyclopedia of genes and genomes pathway and Gene Ontology enrichment analysis
To determine the functions of respective DEGs, we initially performed a Gene Ontology (GO) enrichment analysis (Liu, Liu & Rajapakse, 2018) via the Database for Annotation, Visualization and Integrated Discovery (DAVID) database (Version 6.8, https://david.ncifcrf.gov/) (Huang, Sherman & Lempicki, 2009). The enrichment of different pathways was mapped using the Kyoto encyclopedia of genes and genomes (KEGG) pathway analysis (Kanehisa et al., 2017). An FDR < 0.05 was set as an inclusion criterion in both GO and KEGG enrichment analyses. In order to identify putative roles of distinct gene modules in OSCC, we performed the GO enrichment analysis for DEGs that composed each selected module cluster. The “cluster Profiler” package in R/Bioconductor software (Yu et al., 2012) was utilized for this enrichment analysis. The FDR < 0.05 was set again as an inclusion criterion in this GO analysis.
The included patients’ clinical and pathological data were shown in Table S3. The tissue microarray chips were dewaxed and dehydrated, and then incubated overnight at 4 °C with monoclonal rabbit anti-human SPRR3 (Abcam, ab218131) after epitope retrieval, H2O2 treatment and non-specific antigens blocking, after signals detection by using DAB staining kit (Vector Laboratories, Burlingame, CA, USA), we incubated chips with the secondary antibody. IHC analysis of the chips were performed under optical microscopes of 100× and 400×.
Gene-set enrichment analysis
Gene-set enrichment analysis (GSEA) is a computer-based tool used for microarray data analysis and annotation, which is based on biological knowledge (Subramanian et al., 2005). The GSEA enrichment score (ES) reflects the degree to which a particular set is over-represented at the extremes (top or bottom) of the entire ranking list. ES is normalized for each geneset to account for the overall size of the set, yielding a normalized enrichment score (NES). In this study, the data of OSCC was respectively obtained from our own RNA-seq microarray data and TCGA database. RNA-seq samples were grouped as tumoral and normal tissues. Samples from TCGA database were grouped according to the quantile of SPRR3 expression levels. Thereafter, samples were sorted from low to high in terms of expression level, where the first quantile was defined as the “low expression group” while the last quantile was defined as the “high expression group”. GSEA annotated the GO, KEGG and Hallmark enrichment results of these genes with an enrichment score. The inclusion criteria of the enriched results were established as |NES| > 1, P-value < 0.05 and FDR < 0.25 (Subramanian et al., 2005).
To explore the connection between factors and overall survival (OS) of the OSCC patients, univariate Cox regression model was utilized in specimens of GSE42743. Subsequently, we used significant factors to conduct the Least Absolute Shrinkage and Selection Operator (LASSO) method. LASSO regression contributes to variable selection and regularization, at the same time that it fits the generalized linear model. Using this statistical approach, we could also avoid over-fitting at some extent. Because the complexity of LASSO regression was controlled by the coefficient (λ), a less variable model was obtained by performing a penalty ratio, according to their size. A relatively small number of highly correlated indicators was then obtained and, subsequently, a plot was designed according to the partial possibility deviance vs log(λ). The positive factors of the LASSO regression analysis were incorporated into the multivariant Cox regression model, which confirmed the feasibility of these factors as independent predictors of OSCC. We later verified the effect(s) in the TCGA database using this same methodology.
The K–M survival analysis was performed to examine the relationship between the survival time of OSCC patients and the expression level of respective DEGs. OSCC samples from GSE42743 and TCGA were divided into two groups, according to the quartile of expression for each selected gene. As previously, samples were sorted from low to high in terms of expression level for each target gene. First and last quantiles were related to the low and high expression groups, respectively. We calculated and then indicated the hazard ratio (HR) with 95% confidence intervals and log-rank P-value along the plot.
To evaluate the discriminatory accuracy of a selected gene between two groups, its expression level was included in the Receiver Operating Characteristic (ROC) curve analysis. The ROC curve is an essential tool for diagnostic detection or biomarker prediction in a certain disease (Matthews, Ranson & Kelso, 2011). Plots of sensitivity (true-positive fraction) versus 1-specificity (false-positive fraction) were constructed. The value for the area under ROC curve (AUC) corresponds to the ability of one gene to distinguish tumor tissues from adjacent ones. χ2 test was utilized to verify the correlation between gene expression levels and clinicopathologic factors. In this context, samples were also grouped according to the quantile of expression level of a respective target gene. All statistical analyses were conducted by SPSS 25.0 software (SPSS Inc., Chicago, IL, USA). The correlation coefficient between two genes was calculated by Pearson correlation analysis. For expression data comparison, student’s t test (two-tailed) and ANOVA were performed using GraphPad Prism 8 software.
DEGs search and analysis
According to the sequencing data here available, a total of 229 DEGs were identified in OSCC samples, from which 85 genes were up-regulated while 144 were down-regulated. The most common up- (log2FC > 1, FDR < 0.05) and down-regulated (log2FC < −1, FDR < 0.05) genes are listed in Table S4. Both volcano plot and heatmap were generated to show the discrepancy among the expressed genes when comparing OSCC and adjacent normal tissues (Figs. 2A and 2B).
The enrichment results of the DEGs
Based on the GO enrichment analysis, the commonly enriched categories were (i) keratinization, (ii) keratinocyte differentiation, (iii) cell adhesion, (iv) migration and (v) proliferation. In addition, DEGs were enriched in alcohol and drug metabolic processes (Fig. 2C), which have been closely related to both tumorigenesis and development. The most stringent results of GO enrichment analysis, including BP (biological process), MF (molecular function) and CC (cellular component) terms, (the top 100 items of each terms with the lowest FDR) are presented in Table S5.
Kyoto encyclopedia of genes and genomes pathway analysis showed that the majority of DEGs were enriched in (i) drug metabolism, (ii) retinol metabolism, (iii) chemical carcinogenesis and (iv) leukocyte transendothelial migration pathways (Fig. 2D). Detailed results of KEGG enrichment analysis are presented in Table S6.
In order to validate the DEG-related functions and pathways on an unbiased basis, GSEA analysis was performed using previously generated RNA-seq data. DEGs were mostly enriched in (i) epithelial mesenchymal transition (EMT), (ii) fatty acid metabolism, (iii) K-Ras signaling down and (iv) TGF-β signaling. In regard to the GO-BP terms, GSEA showed enrichment of DEGs in (i) DNA replication, (ii) epidermal cell differentiation and (iii) positive regulation of cell activation. In the context of KEGG enrichment, DEGs were largely present in (i) drug metabolic cytochrome P450, (ii) extracellular matrix receptor interaction, (iii) focal adhesion and (iv) tight junction (Figs. 2E–2G). Above all, DEGs were more concentrated in biological processes including (i) epidermal cell differentiation, (ii) metabolic process and (iii) cell adhesion, reiterating a potential correlation with tumorigenesis and cancer development.
Identification of module clusters via PPI network
A global PPI network was designed after establishing the inter-relationships among the annotated DEGs (STRING database search). As a result, 166 nodes with a combined score > 0.4 were obtained (Fig. S1). MCODE plugin was applied for module analysis, where two modules were chosen for further analysis according to the inter-connectivity and density between nodes in PPI networks. Module clusters 1 and 2 consisted, respectively, of eight and ten nodes (Figs. 3A and 3B). All DEGs present in these two particular modules are shown in Table 1.
|MCODE cluster||MCODE score||Gene ID|
Genes from selected module clusters were included in the GO enrichment analysis to properly understand their functions and gene product attributes. In regard to BP, the DEGs from module cluster 1 were correlated with (i) keratinocyte differentiation, (ii) epidermal cell differentiation and (iii) epidermis development. In terms of CC, the DEGs were enriched in cornified envelope and desmosomes, and for MF, the DEGs were enriched in protein binding, bridging and molecular adaptor activity (Figs. 3C–3E). The DEGs in module cluster 2 were mainly correlated with collagen metabolic process and extracellular matrix organization (Figs. 3F–3E). Interestingly, most of the genes here involved (for instance, Matrix Metallopeptidase (MMP) family members) have been associated with tumor invasion and metastasis (Lu et al., 2017; Zou et al., 2019). The MCODE score of module cluster 1 was higher, and the hub genes in this cluster were barely reported in OSCC studies. Hence, the module cluster 1 was preliminarily chosen for subsequent evaluation.
Validation of DEGs in independent OSCC datasets
According to the TCGA database, OSCC accounts for the majority of HNSCC cases. The RNA expression level of DEGs in module cluster 1 was preliminarily analyzed by GEPIA in the HNSCC dataset. Among these DEGs, CSTA (Cystatin A), EVPL (Envoplakin), PPL (Periplakin), SPRR3, TGM1 (Transglutaminase 1) were significantly down-regulated in tumor tissues when compared with their normal counterparts (P < 0.05) (Figs. 4A–4H). To further demonstrate that these five genes are differentially expressed in OSCC, we utilized the GSE30784 and GSE3524 datasets to confirm potential disease-related DEGs. Consistently, these five genes were significantly down regulated in tumor tissues (Figs. 4I–4M for GSE30784; Figs. S2A–S2F for GSE3524).
Identification of DEG-related prognostic value
To properly select the predictive factors associated with the overall survival of OSCC patients, the five above-mentioned genes and clinical disease characteristics were introduced into univariate Cox regression model using the GSE42743 dataset. The most statistically significant DEGs (SPRR3, PPL, TGM1) and N staging were concomitantly examined by LASSO method. As a result, only SPRR3 and N staging were confidently sorted out (Figs. 5A and 5B). To further determine whether low SPRR3 expression could be an independent predictor of OSCC prognosis, a multivariate Cox regression model was executed. As predicted, multivariate Cox regression analysis showed that SPRR3 expression levels (HR = 0.865, 95% CI [0.754–0.992], P-value = 0.037) act as an independent prognostic factor for the OS of OSCC patients (Table 2).
|Characteristics||Univariate analysis||LASSO||Multivariate analysis|
|HR (95% CI)||P-value||HR (95% CI)||P-value|
|0.801 [0.662–0.969]||0.032||(+)||0.865 [0.754–0.992]||0.037|
|1.970 [1.370–2.832]||<0.001||(+)||1.750 [1.202–2.548]||0.003|
T, primary tumor; N, regional lymph nodes; (+), selected by LASSO regression; (−), abandoned by LASSO regression. The bold characters refers to the P-values of indicators that were statistically significant in the regression models.
The SPRR protein family has been reported to be functional in a variety of tumors (Specht et al., 2013). Still, the expression and function of SPRR3 in OSCC remain largely unclear. Thus, we have presently proposed SPRR3 as an OSCC-related gene and attempted to explore its role in tumorigenesis and development of this condition. Therefore, Kaplan–Meier (K–M) survival analysis was carried out using the GSE42743 dataset to verify whether SPRR3 could act as a prognostic marker (log rank P-value < 0.05; Fig. 5C).
Statistical platforms including univariate Cox regression model, LASSO method and Cox multivariate regression model were also approached to validate the SPRR3’s prognostic value, based on the TCGA database. According to the univariate Cox analysis, factors including low SPRR3 expression, complete tumor remission, lymphovascular and perineural invasion, as well as N staging were significantly associated with a poor prognosis. SPRR3 was also screened out by the LASSO method (Figs. 5D and 5E). Multivariate Cox proportional regression analysis revealed that SPRR3 expression levels (HR = 0.544, 95% CI [0.329–0.899], P-value = 0.017) was an independent prognostic factor for the OS of OSCC patients (Table 3). The predictive value of SPRR3 towards OSCC prognosis was confirmed by K–M survival analysis based on TCGA database (log rank P-value < 0.05; Fig. 5F).
|Characteristics||Univariate analysis||LASSO||Multivariate analysis|
|HR (95% CI)||P-value||HR (95% CI)||P-value|
|0.595 [0.371–0.955]||0.032||(+)||0.544 [0.329–0.899]||0.017|
Moderate and Poor
|Primary therapy outcome
|0.282 [0.176–0.453]||<0.001||(+)||0.242 [0.147–0.397]||<0.001|
|2.109 [1.266–3.512]||0.004||(+)||1.690 [0.991–2.881]||0.054|
|1.727 [1.081–2.759]||0.022||(+)||1.475 [0.857–2.538]||0.161|
|1.833 [1.131–2.971]||0.014||(+)||1.250 [0.711–2.194]||0.438|
T, primary tumor; N, regional lymph nodes; (+), selected by LASSO regression. The bold characters refers to the P-values of indicators that were statistically significant in the regression models.
Low SPRR3 expression is correlated with carcinogenesis and progression of OSCC
According to the OSCC patient data in the TCGA database, SPRR3 expression levels were significantly lower than those in oral normal tissues (P-value < 0.0001, Fig. 6A). ROC analysis was then performed to confirm the statistical relevance of these results. For this, AUC was calculated to identify the diagnostic specificity and sensitivity of SPRR3 expression. Based on the GSE30784 dataset and the TCGA database, down regulated SPRR3 levels yielded an AUC of 0.920 and 0.731, respectively (Figs. 6B and 6C). In this regard, SPRR3 had positive diagnostic performance for OSCC patients’ discrimination. Thus, SPRR3 could act as a potential diagnostic indicator for OSCC.
Student’s t-test (two-tailed) and one-way ANOVA were thus performed and showed that low SPRR3 expression was in the subgroups of patients with high alcohol consumption, poor differentiation, non N0 staging, positive lymphovascular invasion, and positive perineural invasion (Figs. 6D–6H). The χ2 test was also utilized to figure out the correlation between SPRR3 expression of and the clinicopathological characteristics of OSCC. Consistent with previous tests, SPRR3 expression was significantly correlated with alcohol consumption (P-value = 0.011), histologic grade (P-value < 0.001), N staging (P-value = 0.005), lymphovascular invasion (P-value = 0.026) and perineural invasion (P-value = 0.016) (Table 4). These results reiterated that SPRR3 may act as a factor related to tumor progression and/or metastasis.
|Age (years)||<60 yr||31||35||0.515|
|HPV infected||Not infected||70||68||0.575|
|Primary therapy outcome||Progressive disease||16||51||0.153|
T, primary tumor; N, regional lymph nodes. The bold characters refers to the P-values of indicators that were statistically significant in χ2 test.
Differential detection of SPRR3 protein in OSCC
Immunohistochemical staining suggests that the SPRR3 protein levels in OSCC tissues vary according to the histologic grade of the cancer (Figs. 6I–6L). Specifically, SPRR3 levels were lower in tumor tissues from patients at higher histologic grade when compared to those at lower grade.
SPRR3 is enriched in various metabolic pathways and negatively associated with OSCC malignant progression
To probe SPRR3-related pathways in an unbiased manner, GSEA analysis was conducted using OSCC data from TCGA. Based on the hallmark signature, SPRR3 was enriched in (i) apical surface, (ii) fatty acid metabolism, (iii) K-Ras signaling down and (iv) xenobiotic metabolism. GO-BP enrichment showed that SPRR3 was also correlated with (i) alcohol metabolic process, (ii) drug metabolic process, (iii) epidermis development and (iv) fatty acid derivative metabolic process. Lastly, KEGG enrichment assigned SPRR3 into (i) drug metabolism cytochrome P450, (ii) glycolysis/gluconeogenesis, (iii) metabolism of xenobiotic by cytochrome P450, (iv) retinol metabolism and (v) VEGF signaling pathway (Figs. 7A–7C).
To further verify the role of SPRR3 in OSCC carcinogenesis and development, we analyzed the relationship between SPRR3 and other genes from the two module clusters previously screened, using the GSE30784 and TCGA-OSCC datasets. Interestingly, SPRR3 was positively correlated with genes from the module cluster 1 (CSTA, DSG1, EVPL, IVL, PPL, TGM1). CSTA has been reported as a potential OSCC biomarker (Hsiao et al., 2017). In contrast, SPRR3 was negatively correlated with a number of DEGs from the module cluster 2, including PLAU (Plasminogen Activator, Urokinase), COL5A3 (Collagen Type V Alpha 3 Chain), SERPINH1 (Serpin Family H Member 1), MMP 13, MMP7, COL4A6 (Collagen Type IV Alpha 6 Chain), and COL3A1 (Collagen Type III Alpha 1 Chain). Coincidentally, PLAU has been reported as an important factor involved in the regulation of cell migration, proliferation and invasion in OSCC (Zou et al., 2019), while MMP7 appears to have an oncogenic role in OSCC invasion and metastasis (Dasgupta et al., 2012).
By accessing the GSE30784 and TCGA-OSCC datasets, we also examined major SPRR3 co-expressing genes, including potential markers, which were closely associated with metabolic processes and cancer progression. As a result, we observed that SPRR3 was positively correlated with some metabolic process markers including ADH1C (Alcohol Dehydrogenase 1C), ADH7 (Alcohol Dehydrogenase 7), ALDH1A1 (Aldehyde Dehydrogenase 1 Family Member A1), ALDH1A3 (Aldehyde Dehydrogenase 1 Family Member A3), CYP2E1 (Cytochrome P450 Family 2 Subfamily E Member 1). Among the EMT- and metastasis-related markers, we have also verified that SPRR3 was positively correlated with CDH1 (Cadherin 1) and TJP1 (Tight Junction Protein 1). Moreover, SPRR3 was negatively associated with CDH2 (Cadherin 2), FN1 (Fibronectin 1), MMP1, MMP9, SNAI1 (Snail Family Transcriptional Repressor 1), SNAI2 (Snail Family Transcriptional Repressor 2), TGFBI (Transforming Growth Factor Beta Induced) and VIM (Vimentin) (Fig. 7D for GSE30784, Fig. 7E for TCGA). These results indicated that SPRR3 could be involved in a series of metabolic processes, and the downregulation of SPRR3 may be potentially linked to the malignant progression of OSCC.
Oral squamous cell carcinoma, a subset of HNSCC, is one of the most common malignancies from the oral cavity worldwide. It is an aggressive tumor that bears a poor prognosis. In fact, OSCC patients at advanced stages, even after intensive treatment, may develop relapses and distant metastases. Therefore, early diagnosis and intervention remain critical for the management of OSCC. In this regard, measuring tumor markers has gradually become an important method for detecting the tumor status and monitoring therapeutic effects in the clinic.
The development and increasing use of tumor biomarkers has become a pivotal part of cancer diagnosis. Likewise, seeking novel tumor-related markers has provided new trains of thought in OSCC treatment (Csosz et al., 2017). Potential biomarkers with high specificity and sensitivity, identified in datasets by distinct Bioinformatic methods, have been continuously explored (Sun et al., 2015). Previous studies have determined a number of target genes by filtering DEGs in datasets (Yang et al., 2017; Zhao et al., 2018). Still, given the differences in specimens and database platforms, potential markers have failed to bear a significant clinical application value, so their verification process have been far from clinical practice. Hence, we expect that DEGs obtained from biopsy-based sequencing, combined with the validation by independent public datasets and tissues analyses, will be more beneficial for the clinical diagnosis and prognosis of the disease. In this study, we have focused on the identification of potential biomarkers which may possess diagnostic and prognostic significance.
Current RNA-seq data has revealed 229 OSCC-related DEGs that were particularly enriched in processes associated with tumorigenesis and malignant progression. By further screening gene modules in a PPI network, a module cluster containing eight down regulated hub genes in OSCC (cluster 1) was selected as a target. The core gene in this module was retrieved by acquiring OS information from different datasets. Univariate Cox regression, LASSO regression and multivariant Cox regression analyses were successively conducted in order to screen the key genes that could independently predict OSCC prognosis, resulting in the identification of SPRR3 as a central hub gene. Thereafter, we conducted K–M survival analysis and further established ROC curves. As a predictive marker, SPRR3 has considerable clinical significance in OSCC diagnosis and prognosis. Based on clinical data, a low expression of SPRR3 was correlated with patients’ alcohol consumption, histologic grade, lymph node metastasis, lymphatic vascular invasion, and perineural metastasis. GSEA and Pearson correlation analyses have been also conducted for gene annotation. Putting together, we concluded that SPRR3 is potentially associated with epidermis development, metabolic processes and tumor progression.
SPRRs are located in an evolutionarily conserved genetic cluster, designated as the epidermal differentiation complex (EDC) in the chromosome 1q21. Epidermis development and epithelial cell differentiation are some of the basic functions of SPRR proteins (Cabral et al., 2001; Candi et al., 2000). In addition, SPRRs family proteins have been reported to be active during the progression of many cancers. In OSCC, over-expressed SPRR2A (Small Proline Rich Protein 2A) can impair distant metastasis of tongue squamous cell carcinoma cells in vivo (Fang et al., 2016). Down-regulation of SPRR3 expression has been reported in esophageal squamous cell carcinoma, esophageal adenocarcinoma, anal and gastric cancers, while it has been strongly upregulated in colorectal and breast cancers (Cho et al., 2010; Kim et al., 2012). Nevertheless, the role of SPRR3 in oral tissues has never been explored. According to χ2 tests, here we observed that low SPRR3 expression is associated with OSCC progression, lymphovascular invasion and perineural invasion. GSEA enrichment analysis showed that SPRR3 is also correlated with K-Ras and VEGF signaling, which are classic pathways in OSCC carcinogenesis and progression (Caulin et al., 2004; Siriwardena et al., 2018).
TJP1 and CDH1 are classic inhibitors of tumor invasion and metastasis (Li et al., 2016). The positive correlation between SPRR3 and these two target genes indicate that SPRR3 might be essential for cell adhesion. Furthermore, SPRR3 was negatively correlated with some mesenchymal markers, such as CDH2, TGFBI, FN1, VIM, MMPs, and transcriptional factors such as SNAI1, SNAI2 (Li et al., 2016). Therefore, SPRR3 appears to work as a modulating factor involved in the maintenance of oral epithelial tissues and also in EMT inhibition. These particular roles may explain why patients under-expressing SPRR3 are expected to have a relatively poor prognosis. However, the SPRR3’s role in suppressing EMT, as well as its mechanism(s) of action in OSCC, require more in-depth investigation.
Gene-set enrichment analysis enrichments analysis indicates that a series of metabolic processes are impacted by SPRR3. This observation is consistent with other results we acquired by χ2 test, where alcohol consumption is correlated with SPRR3 expression. Particularly, alcohol consumption is an important risk factor for OSCC, and the individual’s ability to metabolize alcoholic beverages somehow reflects on its ability to control cancer. Accordingly, ethanol has been classified as a carcinogenic substance in humans (Zygogianni et al., 2011). Additionally, SPRR3 may be involved in the metabolism of other substances, including fatty acid, drugs, and xenobiotics. It has been reported that ADH, ALDH and CYP2E1 are indispensable in some metabolic processes (Jin et al., 2013). Thus. the relationship between SPRR3 and ADH, ALDH, CYP2E1 can confirm, to some extent, the role of SPPR3 in metabolic processes.
Altogether, our current results reflect the close relationship between SPRR3 and OSCC carcinogenesis and development. In addition, this work provides novel insights for follow-up studies, which should better dissect the impact of SPRR3 in OSCC, at the molecular and cellular level.
In our research, the expression of SPRR3 was detected to be dysregulated in OSCC by means of RNA-seq, and further was verified by Bioinformatics methods. SPRR3 could be a potential bio-marker for identifying OSCC, and its under-expression could predict the prognosis of patients with OSCC. Above all, SPRR3 could be a qualified diagnostic/prognostic biomarker for OSCC.
Protein-protein interaction network.
PPI network of the DEGs identifed from high-throughput microarrays data was constructed. A total of 166 nodes with a combined score > 0.4 were included in the network
Transcriptional level of DEGs in OSCC from GSE3524 dataset.
(A) The expression level of CSTA was significantly lower in tumor group than that in normal group by student’s t test (two tailed) in GSE3524 dataset, P-value = 0.0275. (B) The expression of PPL was significantly lower in tumor group by student’s t test (two tailed), P-value = 0.0200. (C) The expression of EVPL was significantly lower in tumor group by student’s t test (two tailed), P-value = 0.0035. (D) The expression of SPRR3 was significantly lower in tumor group by student’s t test (two tailed), P-value = 0.0011. (E) The expression of TGM1 was significantly lower in tumor group by student’s t test (two tailed), P-value < 0.0001.