Identification of SPRR3 as a novel diagnostic/prognostic biomarker for oral squamous cell carcinoma via RNA sequencing and bioinformatic analyses

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, log₂ 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 log₂ 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.

IHC staining

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, H₂O₂ 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).

Statistical analysis

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. χ² 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- (log₂FC > 1, FDR < 0.05) and down-regulated (log₂FC < −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.

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).

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).

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 N₀ staging, positive lymphovascular invasion, and positive perineural invasion (Figs. 6D–6H). The χ² 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.

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.