Characterization of Cancer Stem Cell Characteristics and Development of a Prognostic Stemness Index Cell-Related Signature in Oral Squamous Cell Carcinoma

Transcriptome profiles by RNA sequencing (RNA-seq; FPKM values) of 328 OSCC tissues were obtained from The Cancer Genome Atlas (TCGA; https://tcga-data.nci.nih.gov/tcga/↗) on March 11, 2020. FPKM values were transformed into TPM values. Meanwhile, the matched clinical characteristics including age, gender, grade, stage, T, M, and N as well as survival information were also retrieved from TCGA database (Supplementary Table 1). By the Ensemble database (http://asia.ensembl.org/index.html↗), the Ensemble IDs were converted into gene symbols.

As described by Malta et al., mRNAsi of each OSCC sample was calculated using a one-class logistic regression machine learning algorithm (OCLR) [9]. mRNAsi was expressed utilizing β values that were ranged from 0 (without mRNA expressions) to 1 (complete mRNA expressions). mRNAsi was retrieved through the multiple platform analyses according to the previous study [9]. Patients were stratified into high- and low-mRNAsi groups on the basis of median mRNAsi. The prognostic value of mRNAsi was determined, and overall survival (OS) between subgroups was analyzed via Kaplan-Meier curves as well as log-rank tests using “survival” and “surviminer” packages.

By employing Estimation of STromal and Immune cells in MAlignant Tumours using Expression data (ESTIMATE) algorithm, the fractions of stromal cells as well as immune cells were inferred in OSCC tissues using gene expression signatures. Based on stromal and immune scores, tumor purity was then estimated. Comparisons of immune/stromal scores and tumor purity in high- and low-mRNAsi groups were achieved with Wilcoxon tests.

The CIBERSORT method (http://cibersort.stanford.edu/↗) was utilized to quantify 22 tumor-infiltrating immune cell fractions in each OSCC tissue based on gene expression profiling [12]. The immune cells included naive B cell, memory B cell, plasma cell, CD8+ T cell, naive CD4+ T cell, CD4+ resting memory T cell, CD4+ memory activated T cell, T cell follicular helper, T cell regulator, T cell gamma delta, resting/activated natural killer cell, monocyte, M0/M1/M2 macrophage, resting/activated dendritic cell, resting mast cell, activated mast cell, eosinophil, and neutrophil. Significant results (p < 0.05) were screened for further analyses.

With “limma” algorithm, DEGs were screened between high and low-mRNAsi samples [13]. Genes with ∣fold change | >1.5 as well as false discovery rate (FDR) < 0.05 were abnormally expressed. These DEGs were visualized into volcano plots.

Using “clusterProfiler” package, functional annotation analyses of mRNAsi-related DEGs were presented for investigating and visualizing biological functions and pathways involving above genes, including Gene Ontology (GO) function annotations as well as Kyoto Encyclopedia of Genes and Genomes (KEGG) [14]. Terms with FDR < 0.05 were considered significantly enriched.

Through univariate Cox regression analysis, DEGs that were closely related to OSCC survival were screened. These prognosis-related key DEGs were optimized through least absolute shrinkage and selection operator (LASSO) analyses with “glmnet” package [15]. The variable selection and the shrinkage of prognosis-related candidate DEGs were then carried out. Afterwards, a risk score model was generated by multivariate Cox regression analysis based on the coefficient and expression of each gene. According to the median value, OSCC patients were clustered into two subgroups. Survival status was observed in two groups. Moreover, OS analyses in patients with high and low risk were achieved by Kaplan-Meier curves. Differences in OS were estimated with log-rank tests. The expression of the mRNAsi-associated genes was visualized into a heat map via pheatmap package. Differences in clinicopathological features (age, sex, grade, and stage) were estimated in high- and low-risk subgroups. The predictive accuracy of the signature was assessed via the area under the curve (AUC) of receiver operator characteristic curve (ROC). Uni- and multivariate Cox analyses were achieved for estimating the predictive independency of this signature. The expression of the mRNAsi-associated genes was detected in 519 OSCC and 44 normal tissues from TCGA dataset. Kaplan-Meier survival curves were conducted for evaluating the prognostic significance of the mRNAsi-associated genes among OSCC patients in TCGA dataset.

The potential molecular mechanisms underlying the mRNAsi-associated signature were carried out via the GSEA method [16]. Enriched KEGG pathways were explored for high- and low-risk samples. Pathways with FDR < 0.05 were significantly enriched.

Factors that were independently predictive of survival were incorporated to construct a prognosis-associated nomogram for investigating the probability of 1-, 3-, and 5-year OS in OSCC patients with glmnet package [15]. The predictive accuracy of the nomogram was estimated by ROC curves. The calibration curves were plotted for observing the nomogram-predicted and observed survival probabilities. Meanwhile, decision curve analyses (DCA) were presented for calculating the clinical net benefit of each factor in comparison to all or none factors [17]. The best model was the one with the highest net benefit as calculated.

Three pairs of fresh OSCC and adjacent normal tissues were collected from OSCC patients who received surgery in the Department of Stomatology, Chinese PLA General Hospital. None of the patients experienced chemoradiotherapy prior to surgery. This study protocol gained the approval of the Ethics Committee of Chinese PLA General Hospital, with written informed consent acquired from each patient (2020-039).

OSCC cells Cal-27 and HSC-3 (Shanghai Cell Bank of Chinese Academy of Sciences (China)) were maintained in DMEM (Gibco, USA) containing 10% fetal bovine serum (FBS). These cells were cultivated in a constant temperature box containing 5% CO₂ and 37°C. They were cultivated onto 6-well plates (5 × 10⁵ cells/well). SiRNA targeting H2AFZ (si-H2AFZ; 5′-CACCGAGACGCTCGATGACTCCGC-5′ (sense), 5′-AAACGCGGAGTCATCGAGCGTCTC-3′ (antisense); RiboBio, China) and its negative control (si-NC) as well as pcDNA3.1-H2AFZ (RiboBio, China) and empty vector were transfected into cells via Lipofectamine™ 2000 transfection kit (Life Technologies, USA). After 48 h, H2AFZ expression was measured by real-time quantitative polymerase-chain reaction (RT-qPCR) and western blot.

Total RNA was extracted from transfected Cal-27 and HSC-3 cells according to the Trizol kit (TIANGEN, USA). The concentration and purity of total RNA were measured. Through GoScript reverse transcription system kit (Promega, USA), total RNA was used to obtain cDNA by reverse transcription reaction. The primer sequences were as follows: H2AFZ: 5′-GGCGGTAAGGCTGGAAAGG-3′ (forward), 5′-TGTCGATGAATACGGCCCAC-3′ (reverse); GAPDH: 5′-GGAGCGAGATCCCTCCAAAAT-3′ (forward), 5′-GGCTGTTGTCATACTTCTCATGG-3′ (reverse). The qPCR reaction was then carried out. The PCR reaction conditions contained predenaturation at 95°C lasting 10 min, denaturation at 95°C lasting 15 s, annealing at 60°C lasting 20 s, and extension at 72°C lasting 20 s, a total of 40 cycles. Afterwards, H2AFZ as well as GAPDH expressions were automatically read and generated in the PCR instrument. The relative expressions of H2AFZ were determined utilizing 2^–ΔΔCt methods.

Cal-27 and HSC-3 cells were seeded onto 96-well plates (5000 cells/well). At the indicated time points (0, 24, 48, and 72 h), 10 μL CCK-8 reagent (Dojindo, Japan) was added to cells and incubated for 1 h. The absorbance at 450 nm was detected for reflecting cell viability using a spectrophotometer.

Cal-27 and HSC-3 cells following transfection were seeded onto 6-well plates (1 × 10⁵ cells/well). After 24 h, the cells were trypsinized without ethylenediaminetetraacetic acid. Then, cells were centrifuged and were treated with 5 μL Annexin V-FITC (Beyotime, China) as well as 10 μL PI at room temperature in the dark for 15 min. BD FACSCalibur Flow cytometry was applied to detect cell apoptosis. Finally, results were evaluated with CellQuest software (BD Biosciences).

Cell invasion was measured with 8 μm Transwells. The membrane was coated with Matrigel. Cal-27 as well as HSC-3 cells following transfection were seeded onto the upper chambers (1 × 10⁵ cells/well). DMEM plus FBS was added into the lower chambers. Following 24 h, the cells from the upper surface of the membrane were scraped and removed using cotton swab. The invaded cells were stained with 0.4% crystal violet. Five fields were counted under an inverted microscope (Olympus, Japan).

Total proteins were extracted from OSCC and normal tissues as well as Cal-27 and HSC-3 cells. Then, BCA kit (Beyotime, China) was utilized for measuring the extracted protein. Protein was then separated with 10% sodium dodecyl sulfate-polyacrylamide gel electrophoresis gels and transferred to PVDF membrane. After being blocked using 5% nonfat milk, membranes were incubated by primary antibodies targeting H2AFZ (1 : 1000; ab214730, Abcam, USA), N-cadherin (1 : 1000; ab98952, Abcam, USA), Vimentin (1 : 1000; ab137321, Abcam, USA), and GAPDH (1 : 1000; ab8245, Abcam, USA) at 4°C overnight, followed by secondary antibodies. The protein bands were developed using enhanced chemiluminescence kit and were quantified by ImageJ software.

R version 3.6.1 (https://www.r-project.org/↗) and GraphPad Prism software version 8.0.1 were applied for statistical analysis. The R packages were retrieved from Bioconductor (http://www.bioconductor.org/↗). Data are displayed as mean ± standard deviation. Wilcoxon test, Student's t-test, and one-way analysis of variance were employed for comparing differences between subgroups. p < 0.05 was indicative of statistical significance.

This study quantified CSCs by mRNAsi scores for OSCC patients from TCGA cohort (Supplementary Table 2). According to the median value, these patients were clustered into high- and low-mRNAsi groups. Survival analyses were carried out for investigating the survival differences in subjects with high- and low-mRNAsi scores. In Figure 1(a), more undesirable survival outcomes were found in the high-mRNAsi group compared to the low-mRNAsi group. Through the ESTIMATE method, we evaluated stromal/immune scores as well as tumor purity in OSCC specimens. There were increased stromal score and immune score as well as lowered tumor purity in the high-mRNAsi group than the low-mRNAsi group (Figure 1(b)). Tumor-infiltrating immune cells were inferred utilizing the CIBERSORT method. Increased infiltration levels of CD8+ T cell, activated memory CD4+ T cell, and resting NK cell were found in high-mRNAsi specimens (Figure 1(c)). Meanwhile, there were higher infiltration levels of resting memory CD4+ T cell in low-mRNAsi specimens.

To identify mRNAsi-associated genes of OSCC, genes with ∣fold change | >1.5 and FDR < 0.05 were screened (Supplementary Table 3). In Figure 2(a), 833 mRNAs displayed upregulation and 197 mRNAs displayed downregulation in high-mRNAsi specimens than low-mRNAsi specimens. Furthermore, these mRNAsi-associated genes were distinctly enriched in tumor-related pathways such as ECM-receptor interaction, proteoglycans in cancer, and PI3K-Akt pathway (Figure 2(b)). GO annotation results demonstrated that these genes were distinctly associated with extracellular matrix organization (Figure 2(c)).

Through univariate Cox regression analyses, 14 mRNAsi-associated genes were distinctly correlated to prognosis of OSCC patients (Table 1). With the LASSO method, 11 optimal candidate biomarkers were selected for constructing a prognosis risk score model (Figures 3(a) and 3(b)). The risk score of each OSCC patient was determined, as follows: H2AFZ expression∗0.0150584034611393 + KPNA2 expression∗0.171144144263817 + CCDC92 expression∗(−0.08377988547353) + GAS1 expression∗(−0.0247182891302941) + NPM3 expression∗0.0519125455291411 + CCL22 expression∗(−0.0581447729360218) + TSPAN11 expression∗(−0.0616269466917819) + CLEC3B expression∗(−0.033743441760667) + TWIST2 expression∗(−0.127555711866237) + TPSAB1 expression∗(−0.0345953977371839) + IGLV2 − 14 expression∗(−0.0208527855597802). All patients were clustered into two groups based on the median mRNAsi (Figure 3(c)). More death cases were found in the high-mRNAsi group than the low-mRNAsi group (Figure 3(d)). In Figure 3(e), high-mRNAsi patients displayed more unfavorable survival than low-mRNAsi patients. NPM3, H2AFZ, and KPNA2 were upregulated in low-mRNAsi samples while CCDC92, GAS1, CCL22, TSPAN11, CLEC3B, TWIST2, TPSAB1, and IGLV2-14 were upregulated in low-mRNAsi samples (Figure 3(f)).

The associations between the mRNAsi-associated signature and clinicopathological characteristics were evaluated in OSCC patients. In Figure 4(a), this signature was markedly associated with stage of OSCC. The AUC of the signature was 0.700, demonstrating the well-predictive accuracy of this signature in predicting survival (Figure 4(b)). As shown in uni- as well as multivariate Cox analyses, age, stage, and risk scores might be independently predictive of OSCC prognoses (Figures 4(c) and 4(d)).

GSEA was carried out to attempt to explain the intrinsic mechanisms underlying the mRNAsi-associated signature. In Figure 5(a), base excision repair, cell cycle, DNA replication, mismatch repair, nucleotide excision repair, ribosome, and spliceosome were distinctly activated in high-risk OSCC. Meanwhile, B cell receptor, calcium, cell adhesion molecules cam, chemokine-chemokine receptor interaction, ECM receptor interaction, focal adhesion, MAPK, and T cell receptor pathways were markedly activated in low-risk samples (Figure 5(b)).

Independent risk factors (age, stage, and risk score) were incorporated to generate a nomogram for prediction of OSCC patients' 1-, 3-, and 5-year OS probabilities (Figure 6(a)). The AUCs at 1-, 3-, and 5-year OS were 0.686, 0.715, and 0.693, confirming the well-predictive accuracy of the nomogram in predicting survival (Figure 6(b)). As demonstrated by calibration curves, nomogram-predicted 1-, 3-, and 5-year OS were highly similar to investigated survival (Figures 6(c)–6(e)). Moreover, DCA results confirmed that higher net benefit of the nomogram for 1-, 3-, and 5-year OS was investigated compared to other clinical factors (Figures 6(f)–6(h)). Hence, this nomogram exhibited the well-predictive efficacy in OSCC prognosis.

The expression of the 11 mRNAsi-associated genes was observed in OSCC and normal tissues. CCDC92, GAS1, H2AFZ, IGLV2, KPNA2, NPM3, and TWIST2 were markedly upregulated in OSCC compared to normal tissues (Figures 7(a)–7(g)). No significant differences in CCL22, TPSAB1, and TSPAN11 were investigated between OSCC and normal tissues (Figures 7(h)–7(j)). Also, CLEC3B displayed decreased expression in OSCC than normal tissues (Figure 7(k)). Among them, H2AFZ expression has been confirmed to be upregulated in several cancer types such as breast cancer [18] and hepatocellular carcinoma [19]. This study proposed the upregulation of H2AFZ in OSCC for the first time. The expression of H2AFZ was further validated in three paired OSCC and normal tissues through western blot. As expected, high expression of H2AFZ was found in OSCC compared with normal tissues (Figures 7(l) and 7(m)). Prognostic significance of the 11 mRNAsi-associated genes was further analyzed across OSCC patients. Our results demonstrated that patients with high expression of CCDC92, CCL22, CLEC3B, GAS1, IGLV2, TPSAB1, TSPAN11, and TWIST2 exhibited the prominent survival advantage (Figures 8(a)–8(h)). In contrast, high expression of H2AFZ, KPNA2, and NPM3 indicated poorer prognosis compared with their low expression (Figures 8(i)–8(k)).

The biological implications of H2AFZ were observed in two OSCC cells. H2AFZ expression was distinctly suppressed by si-H2AFZ as well as was markedly elevated by pcDNA3.1-H2AFZ in Cal-27 as well as HSC-3 cells (Figures 9(a) and 9(b)). CCK-8 assay was carried out for investigating the cell viability of Cal-27 and HSC-3 cells after transfection. Our data demonstrated that si-H2AFZ markedly suppressed cell viability compared to si-NC for Cal-27 as well as HSC-3 cells (Figures 9(c) and 9(d)). Meanwhile, cell viability of Cal-27 as well as HSC-3 cells was significantly heightened under transfection with pcDNA3.1-H2AFZ than empty vector. By flow cytometry, apoptotic levels of Cal-27 as well as HSC-3 cells were observed. As a result, increased apoptotic levels were found in OSCC cells following transfection with si-H2AFZ in comparison to si-NC (Figures 9(e)–9(h)). Also, pcDNA3.1-H2AFZ transfection markedly lowered apoptotic levels of Cal-27 as well as HSC-3 cells than empty vector.

Transwell assays were performed for observing invasive ability of OSCC cells. In comparison to si-NC, the number of invasive Cal-27 as well as HSC-3 cells with si-H2AFZ transfection was markedly decreased (Figures 10(a)–10(c)). Inversely, pcDNA3.1-H2AFZ transfection significantly enhanced invasive ability of Cal-27 as well as HSC-3 cells in comparison to empty vector. Western blot showed that H2AFZ protein was distinctly decreased by si-H2AFZ as well as was distinctly increased by pcDNA3.1-H2AFZ in Cal-27 cells (Figures 10(d) and 10(e)). Also, EMT-related proteins were detected by western blot. We observed that N-cadherin as well as Vimentin expressions were distinctly lowered in Cal-27 cells transfected with si-H2AFZ (Figures 10(f) and 10(g)). However, their expressions were markedly elevated in Cal-27 cells under transfection with pcDNA3.1-H2AFZ. The similar findings were confirmed in HSC-3 cells (Figures 10(h)–10(j)). We further investigated the morphology changes of Cal-27 as well as HSC-3 cells after si-H2AFZ or pcDNA3.1-H2AFZ transfections. Our results showed that compared with si-NC, Cal-27 and HSC-3 cells transfected with si-H2AFZ were in the form of epithelioid carcinoma and the cells became loose (Figure 10(k)). Moreover, after transfection with pcDNA3.1-H2AFZ, Cal-27 and HSC-3 cells changed from epithelioid carcinoma to spindle shape, and the cells became loose. Thus, H2AFZ upregulation promoted invasion and EMT in OSCC cells.