Immunohistochemical Expression of Five Protein Combinations Revealed as Prognostic Markers in Asian Oral Cancer

Oral squamous cell carcinoma (OSCC) has a high mortality rate (∼50%), and the 5-year overall survival rate is not optimal. Cyto- and histopathological examination of cancer tissues is the main strategy for diagnosis and treatment. In the present study, we aimed to uncover immunohistochemical (IHC) markers for prognosis in Asian OSCC. From the collected 742 synthetic lethal gene pairs (of various cancer types), we first filtered genes relevant to OSCC, performed 29 IHC stains at different cellular portions and combined these IHC stains into 398 distinct pairs. Next, we identified novel IHC prognostic markers in OSCC among Taiwanese population, from the single and paired IHC staining by univariate Cox regression analysis. Increased nuclear expression of RB1 [RB1(N)↑], CDH3(C)↑-STK17A(N)↑ and FLNA(C)↑-KRAS(C)↑were associated with survival, but not independent of tumor stage, where C and N denote cytoplasm and nucleus, respectively. Furthermore, multivariate Cox regression analyses revealed that CSNK1E(C)↓-SHC1(N)↓ (P = 5.9 × 10–5; recommended for clinical use), BRCA1(N)↓-SHC1(N)↓ (P = 0.030), CSNK1E(C)↓-RB1(N)↑ (P = 0.045), [CSNK1E(C)-SHC1(N), FLNA(C)-KRAS(C)] (P = 0.000, rounded to three decimal places) and [BRCA1(N)-SHC1(N), FLNA(C)-KRAS(C)] (P = 0.020) were significant factors of poor prognosis, independent of lymph node metastasis, stage and alcohol consumption. An external dataset from The Cancer Genome Atlas HNSCC cohort confirmed that CDH3↑-STK17A↑ was a significant predictor of poor survival. Our approach identified prognostic markers with components involved in different pathways and revealed IHC marker pairs while neither single IHC was a marker, thus it improved the current state-of-the-art for identification of IHC markers.


INTRODUCTION
Head and neck squamous cell carcinoma (HNSCC) is the sixth most common cancer globally (Bray et al., 2018). Every year, more than 700,000 new cases of HNSCC are diagnosed and 350,000 related deaths are reported worldwide. Oral squamous cell carcinoma (OSCC) is the most common cancer of the head and neck region and has a high mortality rate. However, little improvement has been made in the five-year overall survival rate over the years (Bray et al., 2018). Identifying reliable prognostic factors remains challenging. OSCC is believed to originate from the multistep accumulation of heterogeneous genetic changes in squamous cells. These changes progressively enable transformed cells to proliferate and invade (Oliveira and Ribeiro-Silva, 2011). These accumulated changes may explain why tumors at the same clinical stage and localization often show significant differences in clinical outcome.
The main causes of OSCC in Taiwan and some South Asian countries (Belcher et al., 2014) are the consumption of alcohol, tobacco, and betel nut. This contrasts with human papillomavirus (HPV)-positive oropharyngeal SCC which is associated with HPV infection, with higher proportions in western populations than Asian populations (Gillison et al., 2000).
Unlike other malignancies, the relationships between mutations of genes and clinical morphological characteristics such as tumor grade in OSCC are obscure, which has impeded the development of personalized medicine. Cytopathological and histopathological examination of cancer tissues remains the main diagnostic and treatment strategy for OSCC. Although immunohistochemical (IHC) staining may be limited by small volumes taken from samples, varying expression with selected antibodies, and partial reliance on subjective perception, IHC staining provides morphological information about protein expression, and it is simple and cost-effective. Moreover, the procedures and guidelines (Wolff et al., 2007;Hammond et al., 2010;Dowsett et al., 2011) for IHC staining are well established, and widely used in clinics.
The primary aim is to identify a panel of IHC prognostic markers for Asian OSCC, to enable the selection of patients best suited for intensive adjuvant therapy in clinics. Most of previous results on IHC prognostic markers in OSCC were mainly based on one protein, few on two proteins or on one pathway and are reviewed briefly as follows. IHC of cyclin D1, MDM2, and γ-catenin were shown to be potential prognostic markers in a study of 55 patients with buccal SCC who regularly chewed betel nut (Peng et al., 2011). In 2005, IHC of cyclin D1 and Rb overexpression combined with p16 underexpression (denoted by cyclin D1↑-Rb↑p16↓) (Jayasurya et al., 2005) and Rb↓-p53↑ (Soni et al., 2005) were shown to be associated with poor prognosis in a cohort of 348 and a cohort of 98 Indian patients with OSCC, respectively. Moreover, simultaneous coexpression of p53, cyclin D1, and EGFR was a significant prognostic factor in a cohort of 140 Japanese patients with oral cancer (Shiraki et al., 2005). P-cadherin was reported to be marginally significantly associated with poor survival in a small cohort (Muzio et al., 2005). About a decade later, CK1 nonexpression (Lin et al., 2014) and expression of BRCA1 and γH2AX (Oliveira-Costa et al., 2014) were shown to be associated with poor overall and disease-specific survival, respectively.
We started with a list of 742 synthetic lethal (SL) gene pairs collected from the literature, which consisted of several oncogenes, tumor-suppressor genes, genome stability and other cancer genes with important functions. Two genes are termed SL genes if a single mutation of either is not lethal, but their simultaneous mutation leads to cell death (Chang et al., 2016). The SL interactions of these collected pairs in various cancer types are validated either with human cancer cell lines (Bryant et al., 2005;Farmer et al., 2005) or by genome-wide RNA interference (RNAi) knockdown (Barbie et al., 2009;Luo et al., 2009). The list of SL gene pairs can be accessed at 1 . SL pairs are shown to be correlated to survival of cancer cells (Kaelin, 2005). In general, the more cancer cells killed, the better cancer patients' survival. Thus, we speculated that SL pairs are relevant to prognosis assessment. We hypothesized that IHC (protein) expression is concordant to its gene expression, and we used gene expression data of Asian OSCC to select an initial panel of SL pairs which are relevant to OSCC, from the collected SL pairs (of all types of cancer). Next, we adopted the rule of frequently co-expressed gene pairs along with prior knowledge of OSCC to select ∼20 genes for IHC staining. IHC staining is conducted because protein is more stable than mRNA, ultimately functions in cells, and IHC is usable in the clinic. We also combined single IHC into 398 distinct IHC pairs. To identify prognostic markers, we applied Cox regression analysis to each single IHC (each combined IHC pair) and the overall survival of patients with OSCC. Previously, we applied this approach to colorectal cancer (Chang et al., 2016) and lung adenocarcinoma (Liu et al., 2004), and both studies successfully uncovered IHC prognostic markers independent of tumor stage, in addition to revealing novel IHC marker pairs where neither single IHC was a marker. Our approach starts with the collected SL pairs, which allows components participating in different pathways to be identified as prognostic marker pairs. This improves the current state of the art for IHC marker discovery, which hitherto has mainly relied on one protein, few on two proteins or one pathway (Oliveira-Costa et al., 2014). After the prognostic markers were revealed, we validated them using OSCC data with HPV(−) from The Cancer Genome Atlas (TCGA) HNSCC cohort. A schematic graph of our method is presented in Figure 1.

Study Population
A total of 163 cases of oral cavity cancers were identified in the Kaohsiung Medical University Hospital. Although this sample size was moderate, it was less than only four of the 20 and more previous studies. Furthermore, we conducted a large scale of IHC study and the sample size was sufficiently large for multivariate Cox regression analysis. The inclusion criteria for this study were as follows: (1) age at diagnosis of 20 years or older; (2) tumor histology of squamous cell carcinoma with FIGURE 1 | Schematic graph of the study approach. Microarray gene expression of 57 oral cancerous and 22 noncancerous tissues selected OSCC-relevant gene pairs from 742 verified synthetic lethal pairs. Twenty-one genes were marked for immunohistochemistry staining. Pairwise combinations of the 29 IHCs followed by a log-rank test and Cox regression models revealed single/paired and combined prognostic markers. grade 1 to grade 3; (3) ICD-9 site code specific for the oral cavity; (4) patients underwent surgical interventions,; and (5) disease was diagnosed between 2012 and 2014. The exclusion criteria were: (1) patients who underwent biopsy without surgery; (2) patients with secondary malignancy; (3) tumor histology of carcinoma in situ; and (4) SCC from nasopharynx, oropharynx, hypopharynx, and larynx.

Statistical Analyses, Tissue Arrays and IHC Staining
In the following, all statistical tests were two-sided except where otherwise specified, and all analyses were conducted in R software (R Core Team, 2019).

Preprocessing of Gene Expression Profiles for Oral Cavity Cancerous Versus Non-cancerous Tissues
Gene expression datasets were selected based on the following parameters: cancerous and noncancerous tissues, no treatments, no metastasis, and Affymetrix chips (up to November 2010). The OSCC gene profiles conforming to the aforementioned criteria were downloaded from GEO. Mutated genes associated with oncogenesis may differ among various ethnic groups (Ding et al., 2008). Therefore, we collected gene expression data from patients of Han Chinese origin [tissues from patients in Taiwan, GSE 25099 (Peng et al., 2011)], which was the same ethnicity as that of IHC and clinicopathological data used here. Gene expression profiles of the 57 OSCC and 22 noncancerous tissues in the dataset were quantile-normalized using "expresso" in R, then for a given gene the log ratio of its expression in each cancer tissue versus that of the averaged non-cancerous tissues was computed.

Inference of the Initial Panel of Relevant SL Gene Pairs (Table 2) Using Microarray Gene Expression Data
For each SL gene pair, the fractions of (up, up), (up, down), (down, up), and (down, down) patterns were computed, where the cutoff value for up and downregulation was 1.5-fold. The pattern fractions were computed using the log ratios of the microarray gene expression data for the 57 patients with OSCC (GSE 25099).

Permutation Test and False-Positive Rates of the Fractions of Paired Gene Expression
To evaluate the statistical significance (P value) of the fractions of (up, up) and (down, up) patterns of each gene pair in Table 2, for each fraction we conducted a permutation test to generate its nonparametric distribution. The total rearrangements of the labels of (57) cancer and (22) noncancerous tissues was equal to 79 22 , from which we randomly chose 10,000 rearrangements.
For each rearrangement, we computed the fraction of a pattern to form its distribution, from which we assessed the P value of an observed fraction. Moreover, we applied the q-value (Storey and Tibshirani, 2003) ("q value" in R) to estimate the false discovery rate (FDR) of the significance of the gene pairs in Table 2.

Selection of Genes From the Initial Panel for IHC Staining
We first selected genes whose fractions of the (up, up) and (down, up) patterns were ≥15%, except ≥25% (more stringent) for the (down, up) pattern of KRAS SL pair, because the mutation rate of KRAS was only ∼2% in OC and there were 200 and more KRAS pairs. Next, we applied prior knowledge to (i) select CDH3 (the top-3) from the top three partners of EGFR and the top-1 [STK17A, relevant to OSCC (Pickering et al., 2013)] and top-3 (CDK6, a tumor suppressor gene) partners of KRAS from the KRAS SL pairs satisfying the above fraction cutoffs, and to (ii) include genes whose fractions were on the borderline of 15%; this included FEN1-RAD54B [involved in nonhomologous DNA end joining repair pathway (Storey and Tibshirani, 2003); 14%], RB1 [relevant to OSCC (Liu et al., 2004;Presson et al., 2011); 12%] and MSH2-POLB (Kang et al., 2009;Tiong et al., 2014;Chang et al., 2016).

Tissue Microarray Preparation
Clinicopathological features of 163 OSCC patients were collected (Table 1), and their representative cancer specimens were randomly selected from H&E-stained sections and confirmed by pathologists (Chun-Chieh Wu and Yi-Ting Chen). Three cancerous and one noncancerous tissue cores (diameter 2 mm) were longitudinally cut from each paraffin block. The tissue cores were mounted with fine steel needles in new paraffin blocks to produce tissue microarrays. This study was approved by the Institutional Review Board and Ethics Committee of Kaohsiung Medical University Hospital and the Institutional Review Board of Academia Sinica [Nos. KMUHIRB-E(I)-20170034 and AS-IRB-BM-16075]. The data was analyzed anonymously, and therefore no additional informed consent was required. All methods were performed in accordance with the approved guidelines and regulations and the waiver for the informed consent had been obtained from the approving committee.

Immunohistochemistry Staining
Patients cancer samples were cut into 4-µm-thick sections and deparaffinized in xylene as previously described (Chang et al., 2016). Endogenous peroxidase activity was quenched with 3% (v/v) H 2 O 2 . The sections were boiled in 10 mM citrate buffer for 20 min to revive the antigens. The tissues were incubated with 21 primary antibodies at room temperature for 30 min then rinsed three times with phosphate-buffered saline (PBS) (Supplementary Table 1) according to the manufacturer's protocol. The tissues were then incubated at 25 • C for 30 min with secondary antibodies and a horseradish peroxidase/Fab polymer conjugate [EnVision detection systems peroxidase/DAB, rabbit/mouse (K5007 HRP; DaKo)] then rinsed three times with PBS. Finally, chromogen was developed using 3,3 -diaminobenzidine tetrahydrochloride as the substrate, and counterstained with hematoxylin and viewed under a microscope. Staining intensity in the cancer tissue was independently examined by two pathologists (Chun-Chieh Wu and Yi-Ting Chen). The scoring criteria used here were the same as those of previous studies (Su et al., 2004;Tiong et al., 2014;Chang et al., 2016) (Supplementary Table 2). Stain intensity is graded as negative (0), indeterminate (±), weakly positive (1+), moderately positive (2+), or strongly positive (3+). The criterion is exactly based on the strongest intensity followed by the % expression of the detected protein. Negative (0) indicates no expression of the detected protein, indeterminate means that the staining is weak and its percentage cannot be accurately counted, weakly positive indicates <5% expression of the detected protein, moderately positive is focal expression in 5-20% of the cancer cells, and strongly positive indicates diffuse expression in >20% of the cancer cells. The mean staining intensity of three cancerous tissues was compared with that of noncancerous oral mucosa and a NA denotes missing data; b T and N denote tumor size and lymph node status of "TNM" (AJCC version 7), respectively.
was categorized as either over-or underexpressing, determining the criteria for IHC analyses. The cutoffs for the 29 IHC stains are listed in Table 3.

Log-Rank Test
For each individual and paired IHC staining, a log-rank test of the "high-" and "low-" risk patients was conducted. The highand low-risk groups consisted of patients classified according to the IHC amounts shown in Table 3. If the log-rank test was significant (P < 0.05), which indicated the survival curves of the two groups significantly different, then a Kaplan-Meier survival curve was plotted by the R software.

Univariate Cox Proportional Hazard (PH) Regression Models
In the univariate Cox regression models, the associations between the 29 individual IHC or 398 combined IHC staining pairs and the 10-year overall survival of the OSCC patients were analyzed in the study cohort. The associations between clinical factors such as age (>55 vs. ≤55 years), sex (male vs. female), tumor grade The p-value for the highest fraction four patterns was computed by permutation test with 10,000 repeats, and the false discovery rate was estimated by q-value. Fractions of the four differentially expressed patterns based on the 1.5-fold threshold and filtered from 742 synthetic lethal gene pairs.
(medium and high vs. low), lymph node metastasis (yes vs. no), stage (III, VI vs. I, II), and habits alcohol use (yes/no), betel nut chewing (yes/no), and cigarette smoking (yes/no) with 10-year Taiwanese OSCC overall survival were also assessed.

Multivariate Cox PH Regression Model
When fitting the multivariate Cox regression models, the clinical factor stage significantly associated with overall survival in the univariate Cox regression models was adjusted, because the stage had stronger significance than that of the grade. Likelihood ratio test (LR χ 2 ) and the statistical significance values generated (P values) were used to compare model fit between the uncovered prognostic IHC markers.

Determination of the Cutoff for Differential Expression of the TCGA Data
We first used 1.5-fold as the threshold for differential TCGA OSCC gene expression, but there were too few patients (less than 5) (Vittinghoff and McCulloch, 2007) in the poor/good overall survival subsets to perform univariate Cox regression for most of the six prognostic markers (Table 4A). Thus, the cutoff was relaxed to 1.4-fold, and there were ensure adequate numbers of patients in the poor/good overall survival subsets of two pairs CDH3-STK17A and FLNA-KRAS, respectively, for the univariate Cox regression analysis.

Description of Study Population
As shown in   (Table 2), where up and down denoted upregulation and downregulation with the cutoff 1.5-fold; this less stringent cutoff was set to include important OSCC onco-and tumor suppressor genes not expressed at twofold level, e.g., TP53, EFGR, and CDKN2A, but that were frequently mutated in Asian OSCC (Liu et al., 2004;Presson et al., 2011). Overexpression of tumor suppressor-and genome stability gene pairs associated with DNA repair such as BRCA1 and FEN1 was unexpectedly noted (Table 2). However, this finding was consistent with the dramatic increase in genomic instability and DNA replication caused by mutant oncogenes such as MYC.
Twenty-One Genes Were Selected for IHC Staining We selected 21 genes from Table 2 to conduct IHC staining, and some of them were stained at two cellular portions. Most of the genes were selected according to relatively high fractions of the (up, up) and (down, up) patterns (≥15%) in Table 2. For an extended list of the sorted (up, up) and (down, up) gene pairs, please see 2 . Next, we applied prior knowledge to (i) select CDH3 from the EGFR SL pairs and STK17A and CDK6 from the KRAS SL pairs, which satisfied the above fraction cutoffs, and to (ii) include genes whose fractions were on the borderline of 15%; this included FEN1-RAD54B (Srivastava and Raghavan, 2015) (14%), RB1 (Liu et al., 2004;Presson et al., 2011) (12%) and MSH2-POLB (Kang et al., 2009;Tiong et al., 2014;Chang et al., 2016). Please see the section "Materials and Methods" for details of the selection method. Eight out of these 21 genes were stained at two cellular portions, such as CDH3 and EGFR, the remaining 13 genes were stained at one cellular portion. Table 3  We next explored if the results of IHC stains are suitable for use as prognostic markers. For each of the 29 IHC results in Table 3 and all of the combined IHC pairs, we applied log-rank tests to the 153 Taiwanese patients with OSCC for whom overall survival was recorded. We first observed that the patients with overexpressed RB in nucleus (denoted as RB1(N)↑) had significantly poorer overall survival than patients with underexpressed RB1 (P = 0.027, Figure 3A). Additionally, underexpressed FLNA in cytoplasm [FLNA(C)↓] was also associated with poor clinical outcomes (P = 0.047, Figure 3B).

IHC of Eight Protein Pairs Were Associated With Overall Survival
Furthermore, we combined the 29 IHC stains into all the possible distinct IHC pairs (398 in total), which allowed novel paired IHC markers to be uncovered, excluding those of the same protein stained at different cellular
Of the clinical variables [age, sex, tumor grade, lymph node (LN) metastasis and stage], grade, lymph node metastasis and stage were significantly associated with the patients' overall survival ( Table 4A). The univariate model based on stage (LR 2 x = 13.1) fit better than that based on grade (LR 2 x = 4.8). Therefore, we used stage as the adjustment factor in the multivariate Cox regression models.
Because the high incidence of oral cancer in Asian OSCC is related to alcohol use, betel nut chewing, and cigarette smoking, we investigated whether these habits were associated with overall survival in this population. As shown in Table 4A, only alcohol use [hazard ratio (95% confidence interval) = 2.01 (1.01-3.97); P = 0.045] was a significant predictor of overall survival in these Taiwanese patients with OSCC. Betel nut chewing [hazard ratio (95% confidence interval) = 0.72 (0.38-1.38); P = 0.329] and smoking [hazard ratio (95% confidence interval) = 1.79 (0.64-5.00); P = 0.267] were not significant predictors for the risk of death in these patients. Furthermore, the correlation between alcohol use and each IHC marker was tested (Fisher's exact test) and none was significant at P = 0.05. Similarly, the correlation between "the combined habits" and each IHC marker was tested, and again none was significant at P = 0.05; please see Supplementary Table 3 for details.
We then evaluated the associations of the novel IHC prognostic markers with overall survival after adjusting for alcohol use, age and tumor stage, via multivariate Cox regression analysis. The paired prognostic markers CSNK1E↓-SHC1(N)↓ [hazard ratio (95% confidence interval) = 7.75 (2.85-21.07); P = 5.9 × 10 −5 ], CSNK1E↓-RB1(N)↑ [hazard ratio (95% confidence interval) = 2.16 (1.02-4.58); P = 0.045], and BRCA1(N)↓-SHC1(N)↓ [hazard ratio (95% confidence interval) = 2.87 (1.11-7.42); P = 0.030] were significant predictors of the overall survival of the patients with OSCC. For patients with CSNK1E↓-SHC1(N)↓, CSNK1E↓-RB1(N)↑, and BRCA1(N)↓-SHC1(N)↓, the risk of death was 7.8, 2.2, and 2.9 times higher, respectively, than that for the other patients in this population. However, RB1(N) [hazard ratio (95% confidence interval) = 1.71 (0.89-3.30); P = 0.108] was no longer a significant predictor (Table 4B). After alcohol use, age, and stage were entered into the multivariate Cox regression models along with each of the markers, neither alcohol use nor age was selected by stepwise selection or Akaike information criterion (AIC). Thus, neither appeared in the final models (Table 4B). Following a reviewer's suggestion, we further adjusted the effect of lymph node metastasis and tumor stage in the multivariate Cox regression models, because lymph node density and metastasis were shown to be significant prognosis predictors in OSCC (Zanaruddin et al., 2013;Chang et al., 2018). Excluding the effect of LN metastasis and stage that are used in clinical practice conventionally, the revealed five combined markers are still significant (Supplementary Table 4). This highlights the potential of these markers being targeted for cancer treatments.

Combinations of Significant Markers Were Studied; CSNK1E↓-SHC1(N) Was Suggested for Clinical Practice
We then combined any two of the significant markers in Table 4A and selected eight combinations whose good/poor OS subsets consisted of a sufficient number of (≥5) patients. Note that the (↑,↑) subset of FLNA(C)-KRAS(C) was correlated with a good OS, so we combined its complementary subsets (↓, ↑), (↑, ↓), and (↓, ↓) with the poor OS subsets of the remaining five markers in Table 4A. We fitted multivariate Cox regression 29 TP53(N) ≤0% >0% a The notation (C), (N) and (M) represent cytoplasm, nucleus and membrane, respectively. b Both strength of staining ≥ 1+ and stained area ≥ 10% are required. c Phosphorylated MET was stained. d The stained area > 0% is required.
to these combinations, and found that [CSNK1E-SHC1(N), FLNA(C)-KRAS(C)] (hazard ratio = 8.71; P = 0.000 rounded to three decimal places) and [BRCA1(N)-SHC1(N), FLNA(C)-KRAS(C)] (hazard ratio = 3.14; P = 0.020) were significant prognostic markers (Table 4C). For combinations of three or more significant markers, the (good/poor OS) subsets had too few patients to fit any multivariate Cox regression model. Taking Table 4A-C together, we suggest using the combination CSNK1E↓-SHC1(N)↓, which has the most significant P value (from likelihood-ratio test) among all markers, to identify Asian OSCC patients with worst survival in clinical practice.

External Validation of the Association of CDH3-STK17A With Overall Survival
Ethnicity and geography may play a role in the etiology of cancer. If the newly discovered markers are confirmed by independent  datasets of patients with OSCC from different geographic regions and ethnicities, they may be useful tools in clinical medicine. OSCC with HPV(−) from the TCGA head and neck SCC cohort (henceforth, TCGA) (The Cancer Genome Atlas Network, 2015) more closely resembled Asian OSCC than those with HPV(+). Therefore, we analyzed microarray gene expression data of 160 OSCC cases with HPV (−) to validate the novel prognostic markers in Table 4. We used 1.4-fold as the cutoff for differential expression of TCGA OSCC RNA-seq data (see section "Materials and Methods" for details), such that of all markers in Table 4A, CDH3-STK17A and FLNA-KRAS had a sufficient number of (five or more) (Vittinghoff and McCulloch, 2007) patients in the (good/poor OS) subsets for the univariate Cox regression analysis. Of these, CDH3↑-STK17A↑ was a significant gene predictor of good survival [hazard ratio (95% confidence interval) = 0.55(0.35-0.87); P = 0.011], while FLNA↑-KRAS↑ was not significant (P = 0.117). The former finding was not consistent with ours, which showed that CDH3(C)↑-STK17A(N)↑ was a significant predictor of poorer overall survival in Taiwanese patients with OSCC (Table 4A). This discrepancy may be explained by the different genetic backgrounds in the two populations, since the significant downregulation of the CDH3 gene has been reported in metastatic OSCC (Méndez et al., 2009). The estimated survival curves of CDH3-STK17A are shown in Figure 4. This external validation demonstrates that if IHC or gene expression data are available, CDH3-STK17A can be used to stratify patients with OSCC in the future.

DISCUSSION
Here, we established a cost-effective approach for the identification of prognostic IHC markers of OSCC. This approach is also efficient, as merely 29 IHC stains were performed, but five clinically beneficial prognostic markers were identified through extensive statistical analysis. Our technique rapidly uncovered the prognostic markers without any prerequisite knowledge of the molecular pathways. In contrast, previous studies relied on pathway information (Sadanandam et al., 2013;Kosari et al., 2014) to reveal prognostic markers, such as cellular phenotypes and protein expression levels. Moreover, our approach was able to reveal IHC prognostic markers with components from different pathways. This improves the current state of art, as most of methods to uncover IHC markers to date have been mainly based on one or two proteins (Lin et al., 2014), or one pathway (Oliveira-Costa et al., 2014).
Of the single IHC results, RB1(N) was a predictor of poorer survival in the Taiwanese patients with OSCC, however, it was not independent of tumor stage. This finding was consistent with earlier studies wherein RB1 was a biomarker in HPV(−) head and neck cancers (The Cancer Genome Atlas Network, 2015;Beck et al., 2016). Previous studies showed that expression of Rb increased in the development and/or with disease progression of OSCC (Pavelic et al., 1996;Schoelch et al., 1999;Thomas et al., 2015), and the latter study reported high expession of Rb in patients with combined habits (alcohol use, betel nut chewing and smoking), suggesting Rb pathway altered. However, in our study, the over-expression of Rb was confounded with stage (Tables 4A,B), but not associated with the combined habits (P = 0.19, Fisher's exact test; Supplementary Table 3). The overexpression of RB1(N) in our study may be due to over-expression of cyclin D1 or under-expression of p16 INK4A , as cyclin D1 and p16 INK4A are related to Rb through an autoregulatory loop (Gimenez-Conti et al., 1996;Andl et al., 1998). Although expression of RB1(N) was high in our study, its function was likely inactivated which may be due to regulation of cyclin D1, HPV infection (Gimenez-Conti et al., 1996;Andl et al., 1998), loss of heterozygosity (Maestro et al., 1996;Yokoyama et al., 1996) or Rb hyperphosporylation (Chatterjee et al., 2004), but further studies are required to elucidate this.
Of the 398 IHC pairs, multivariate Cox regression analyses showed that CSNK1E↓-SHC1(N)↓, CSNK1E↓-RB1(N)↑, and BRCA1(N)↓-SHC1(N)↓ were significant predictors of the risk of death in this Taiwanese OSCC population, independent of tumor stage. Of all combinations of two significant markers in Table 4A, , FLNA(C)-KRAS(C)] was the most significant poor prognostic factor. Nevertheless, this marker was less significant than CSNK1E↓-SHC1(N)↓ statistically. Thus, in clinical practice we recommend using CSNK1E↓-SHC1(N)↓ to identify patients with severe and/or advanced Asian OSCC, who should be suggested for alternate or more intense treatment strategies in clinical practice. CK1 could be an oncoprotein or a tumor suppressor (Lin et al., 2014), but phosphorylation of CK1 can stabilize and activate tumor suppressor p53 (Knippschild et al., 2005). SHC1 is a known downstream target of p53, which involves in stress-induced signal transduction pathway (Trinei et al., 2002). Moreover, SHC1 was downregulated by miR-5582-5p, thus led a tumor suppressive activity with GAB1 . In our study, the mean survival rate of OSCC patients with CSNK1E↓-SHC1↓ is 13.8 months compared to 37.8 months of the remaining group. Collectively, CSNK1E-SHC1 might be a tumor suppressor, but this requires further studies for elucidation. As phosphorylation of CK1 can stabilize and activate tumor suppressor p53, moreover, expression of p53 was lower in OSCC lesions than in malignant lesions, and Rb expression was observed in OSCC lesions (Oliveira and Ribeiro-Silva, 2011). Thus, we speculate p53 may indirectly interfere with RB1, after p53 been regulated by phosphorylated CK1, which supports the finding CSNK1E↓-RB1(N)↑ is a poor prognostic marker.
External gene expression data of HPV(−) OSCC from the TCGA cohort (98.7% non-Asian patients) validated that CDH3↑-STK17A↑as a significant predictor of good survival in 160 patients with HPV(−) OSCC, where the cutoff was set at 1.4-fold. Nevertheless, when we set the cutoff at 1.5-fold, this gene pair was no longer significant (P = 0.124). Thus, this gene pair may not be a robust prognosis marker. Our result showed that CDH3(C)↑-STK17A(N)↑ was correlated with poor survival of 163 Taiwanese patients with OSCC. The difference in the aforementioned results may be because overexpression of P-cadherin (coded by CDH3) in membrane was associated with good survival of 67 OSCC patients, however, cytoplasmic expression of P-cadherin was correlated with poor survival (Muzio et al., 2005). Furthermore, high cytoplasmic expression of STK17A was reported to increase tumorigenic potential through inhibition of TGF-beta1-mediated tumor suppressor activity in HNSCC cells (Park et al., 2015).
Some of the prognostic markers our approach revealed are well-known and reported in the literature, thus we performed a comparative analysis as follows. Among all biomarkers, CSNK1E↓-SHC1(N)↓ was the most significant in terms of P < 0.0001 (Table 5). Consistently, the loss of CK1ε expression was shown to be a poor prognostic marker in Taiwanese patients with oral cancer (Lin et al., 2014). Next, overexpression of cyclin D1 and Rb and low expression of p16 was significantly associated with reduced disease-free survival in 348 Indian patients with OSCC (Jayasurya et al., 2005), consistent with our findings that the overexpression of cytoplasmic Rb was a poor prognostic marker in Taiwanese patients. However, Soni et al. (2005) reported that 105 Indian patients with OSCC with loss of Rb expression had poor prognosis. This discrepancy may be explained by a recent finding (Sanidas et al., 2019) that there are many different forms of active Rb, and they have distinct functional properties. Both Shiraki et al. (2005) and our team found that overexpression of p53 (EGFR) was not a significant prognostic marker in OSCC, but the former study revealed that p53-Cyclin D1-EGFR was significantly associated with poor overall survival (P = 0.019). Moreover, BRCA1 overexpression was shown to be associated with reduced overall survival of 150 Brazilian patients with OSCC (Oliveira-Costa et al., 2014), whereas we did not find prognostic significance of BRCA1 underexpression, but BRCA1(N) ↓-SHC1(N)↓ was an independent prognostic marker (P = 0.024; Table 4B). This discrepancy may be due to the different genetic backgrounds of the populations.
In conclusion, our study revealed that the combined evaluation of CSNK1E↓-SHC1(N)↓ in OSCC identified a group of patients with the poorest survival, who should be suggested to undergo alternate or more intense treatment strategies. CK1ε combined with SHC adaptor protein 1 emerged as the most promising IHC prognostic marker in Asian OSCC. Of the 398 combined IHC pairs, genes of ten pairs are known to be SL, out of which only FLNA-KRAS was revealed to be a good OS marker, but not independent of stage. Excluding the effect of tumor stage and LN metastasis (Zanaruddin et al., 2013;Chang et al., 2018) that are used in clinical practice conventionally, the revealed markers of our study are still significant (Supplementary Table 4). This highlights the potential of these markers being targeted for cancer treatments.
Despite that we conducted a large scale of IHC study, the present study is limited by the moderate sample size and no genomic data profiled. Further studies based on larger sample sizes of patients with OSCC and on DNA sequencing data will reveal whether the expression of the uncovered IHC markers are due to their mutations. With ready availability of gene expression and tissue array data and resources to match clinicopathological features in the public and commercial domains, our approach can immediately be applied to other types of cancers. Moreover, additional IHC stain of cyclin D1 will enable us to evaluate the prognostic significance of protein triplets such as cyclin D1-Rb-p16 and p53-cyclin D1-EGFR. This is interesting, as the component of our most significant marker SHC adaptor protein 1-CK1↑ is involved in the EGFR pathway and is SL to both TP53 and EGFR, respectively. Given that the triplet IHC cyclin D1-Rb-p16 is a promising marker, future studies will extend to the prognostic effect of triplets of IHC in OSCC.

DATA AVAILABILITY STATEMENT
Part of the datasets in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found in the articles. However, the clinical data that support the findings of this study are available from Kaohsiung Medical University Hospital with restrictions applying to the availability of these data, which were used under license for the current study, are not publicly available. Data are however available from the authors upon reasonable request and with permission from Kaohsiung Medical University Hospital.

ETHICS STATEMENT
This study involving human participants was approved by the Institutional Review Board and Ethics Committee of Kaohsiung Medical University Hospital (KMUHIRB-E(I)-20170034). The data were analyzed anonymously, and therefore, no informed consent was required. All methods were performed under approved guidelines and regulations.

AUTHOR CONTRIBUTIONS
GS conceived the study. H-CW, C-CW, and Y-TC did IHC staining. C-JC implemented the methods, wrote the algorithm, and performed data analysis. H-CW, J-GC, and GS interpreted the data. H-CW and GS wrote the manuscript; T-CL modified H-CW's writing in an earlier version. GS and T-CL designed and supervised the study. All authors read and approved the final manuscript.