Bioinformatics Analyses Indicate That Cathepsin G (CTSG) is a Potential Immune-Related Biomarker in Oral Squamous Cell Carcinoma (OSCC)

Introduction

About 2.7*10⁵ new oral squamous cell carcinoma (OSCC) cases were reported every year,¹ OSCC is one of the most common cancers all over the world. A majority of OSCC lesions develop from potentially malignant disorders, and are characterized by invasion, distant metastases, high recurrence rate, and undesirable prognosis.² The most common risk factors include smoking, alcohol, and betel nut consumption.³ For the past few years, therapeutic methods of OSCC, including operative treatment, radiotherapy and chemotherapy, have made great progress. However, the postoperative 5-year overall survival rate is only 50%⁴ owing to lymph node metastases, distant metastases, and terminal-staged diagnosis.⁵ Therefore, in order to raise the overall survival rate of OSCC, it is a prerequisite and imperious to search for efficient biomarkers for OSCC.

Recently, many studies have shown that the immune system played a significant role in tumorigenesis.⁶ The immune system could inhibit cancer progression through formatting inflammatory microenvironment and then removing tumor cells and blocking cell proliferation and invasion.^7,8 Furthermore, HPV infection was identified as a risk factor in OSCC,⁹ indicating that a tumor-related inflammatory microenvironment existed in all stages of tumorigenesis. Generally, it was believed that abnormal adaptive and innate immune responses were relevant to the occurrence of tumors by selecting for invasive cloning, inducing immunosuppression and promoting the progression and metastasis of cancer cells.¹⁰ Accumulating evidence indicated that immune evading was related to the initiation and progression of head and neck squamous cell carcinoma (HNSCC).^11,12 In Hepatocellular Carcinoma (HCC),¹³ immunotherapies based on PD-1/PDL-1 showed promising clinical responses. In OSCC, Caldeira et al¹⁴ demonstrated that neutrophil infiltration was different between the carcinoma size T1/T2 group and the T3/T4 group, with higher indexes in the advanced lesions. Shimizu et al¹⁵ identified that CD8+ T-cells in the parenchyma of the invading edge and peripheral stroma were related to tumor recurrence and prognosis. In total, the immune system dysregulation was significantly relevant to the occurrence and progression of OSCC. The immune-related genes, as a significant section of the immune system, might be effective therapeutic targets in OSCC. At present, there were various studies proposing gene expression signatures for predicting overall survival in OSCC patients.^16,17 However, fewer studies paid attention to the effect of whole subsets immune-related genes (IRGs) in OSCC. Given the significant effect of the immune system in the pathogenesis of OSCC, our study intended to reveal the prognostic value of IRGs and screen the effective immune-related biomarkers for OSCC by integrated bioinformatics analysis.

In this study, we downloaded RNA-seq data and relevant clinical information from The Cancer Genome Atlas (TCGA) and extracted IRGs and transcription factors (TFs) expression profiles from the datasets. The present study aimed to analyze the differentially-expressed IRGs and TFs, build a co-expression network between IRGs and TFs, select differentially-expressed IRGs to construct a gene expression profile signature, and identify crucial prognostic immune-related biomarkers in OSCC samples.

Materials and Methods

Data Acquisition and Differentially-Expressed Analysis

IRGs were obtained from the ImmPort database (https://www.immport.org/) and TFs downloaded from the Cistrome Data Browser (http://cistrome.org/). OSCC RNA-seq data were downloaded from TCGA (https://portal.gdc.cancer.gov/) following extracting IRGs and TFs expression profiles. The OSCC samples with incomplete clinical data were excluded. Subsequently, EdgeR in R software was used to screen differentially-expressed IRGs with the standard |log₂ (fold change [FC])|>1.0 and FDR (adjusted P-value)<0.05.

Samples Collection

Forty-eight OSCC cases and their normal controls (matched normal tissues, MNTs) were collected from NanFang Hospital, Southern Medical University. All participants (including healthy controls) have written and provided informed consent to participate. All cases were pathologically confirmed to be squamous cell carcinoma and normal tissues, respectively.

Construction of Co-Expression Network

Differentially-expressed IRGs and TFs were enrolled in the co-expression network. The regulatory network was constructed by the “WGCNA” package in R software with the criteria correlation ≥0.4. And the network was visualized in Cytoscape software.

Establishment of Gene Expression Signature

All differentially expressed IRGs were selected to screen prognostic factors through univariate Cox regression analysis with the cut-off standard P-value<0.05. In order to construct a precise gene expression signature, these prognostic IRGs were selected into multivariate Cox regression and step-wise regression analysis. Finally, the risk level of each OSCC patient was calculated and visualized in R software according to the risk Cox formula and gene expression level.

Correlation Analysis

Cox regression was used to analyze the clinical characteristics, including risk score, clinical stage (stage Ⅰ, Ⅱ, III, Ⅳ), age, tumor grade (grade Ⅰ, Ⅱ, III, Ⅳ), T, N classification in TNM system, and gender. Distant metastasis (M classification) was excluded owing to numerous missing data. Subsequently, correlation analyses between risk level and these clinical characteristics were explored. In addition, the immune cells level was obtained from Tumor IMmune Estimation Resource (TIMER) (https://cistrome.shinyapps.io/timer/), and correlation analysis of immune cells level and risk score was also performed in R software.

Survival Analysis

On the basis of gene expression signature, OSCC samples were classified into low- and high-risk level groups. Kaplan–Meier survival curves, combined with stratification analysis, along with a log-rank P-test were performed to verify its accuracy in the R software “survival” package.

Cell Culture and Transfection

OSCC cell lines including SCC9, SCC15, SCC25, and normal Human Oral Keratinocytes (HOK) were purchased from the Institute of Antibody Engineering, Southern Medical University (Guangzhou, China). And the use of all cell lines was approved by the Nanfang Hospital ethics committee. HOK, SCC15, and SCC25 were cultivated in DMEM (Gibco, Cat#11995500TB) and SCC9 in DMEM/F12 (Gibco, Cat#C11330500BT) with 10% fetal bovine serum (FBS) (ExCell Bio, Inc., Shanghai, China) and incubated at cell incubator. The expression vectors for CTSG were designed and synthesized in Tsingke (Guangzhou, China). Expression vectors were transfected into OSCC cells following the protocol of lipofectamine 2000 (Thermo, Cat# 11668019). Total RNA and protein were extracted after 48–72 hours.

RNA Extraction and Quantitative Real-Time PCR (qRT-PCR)

Total RNAs were extracted in a TRIZOL (Takara, Cat# 9109) manner. Then, equal total RNA was reversed to cDNA on the basis of the manufacturer’s protocol (Vazyme Biotech Co., Ltd, Cat# R212-02). The relative expression levels of target genes in OSCC cases were quantified by RT-qPCR with ChamQ Universal SYBR qPCR Master Mix (Vazyme Biotech Co., Ltd). And GAPDH was used to normalize expression. CTSG Forward primer (5ʹ-3ʹ):GAGTCAGACGGAATCGAAACG; CTSG Reverse primer (5ʹ-3ʹ):CGGAGTGTATCTGTTCCCCTC. GAPDH Forward primer (5ʹ-3ʹ):CGCTGAGTACGTCGTGGAGTC; Reverse primer: (5ʹ-3ʹ):GCTGATGATCTTGAGGCTGTTGTC.

Cell Viability Assay

Approximately 2*10^3 OSCC control cells and CTSG overexpression cells were seeded in 96-well plates, Cell Counting Kit (Vazyme, Cat#A311-02-AA) was used to detect the cell viability at 1, 2, 3, and 4 days, respectively, with Biotek synergy HTX. Ninety microliters serum-free medium and 10 µl CCK8 reagent were mixed and added to each well. Then, these cells were incubated for 2 hours in the cell incubator. Finally, cell viability was reflected by the absorbance at 450 nm.

Clonogenic Assay

A sample of 0.25% trypsase (Solarbio, Cat# T1300) was used to digest the transfected cells and then suspended with complete medium. Then, control and CTSG overexpression cells were both cultured in 6-well plates with the starting density of 2,000 cells per well in a cell incubator for 1–2 weeks. Subsequently, cells were washed with PBS twice and fixed by 4% paraformaldehyde for 10 minutes, following staining for 15 minutes by 0.1% crystal violet staining solution. Finally, an inverted phase-contrast microscope was used to take colony images.

Transwell Assay

The 100l matrigel 1:8 dilution was coated at the bottom of the chamber and then placed in a cell incubator for 1–2 hours. Fifty thousand transfected cells were resuspended with 200l DMEM and seeded in the upper chambers. Then 600 μL 10% complete medium was placed in the bottom wells. The migration assay was without matrigel in the chamber. Subsequently, the transwell device was put in the cell incubator for 2–3 days. Then, these cells in the upper side were removed by cotton swab, and these cells below the chamber were fixed by 4% formaldehyde and then stained with 0.1% crystal violet. Finally, migration and invasion cells were viewed in an inverted microscope.

Western Blotting

Cells and tissue protein were extracted by RIPA lysis. SDS PAGE gels were used to separate protein samples, and then proteins were transferred on PVDF membranes (Pall, Cat# BSP0161). Subsequently, PVDF membranes were sealed with 5% skim milk for 1 hour and then primary antibodies were incubated at 4°C overnight. Goat anti-rabbit secondary antibodies for 1 hour at room temperature, and finally proteins were detected by ECL (YEASEN, Cat# 36208ES76) Prime Western Blotting Detection reagent. CTSG (ZEN BIO, Cat# 513172,1:1,000); E-cadherin (Proteintech, Cat# 20648-1-AP,1:3,000); N-cadherin (Proteintech, Cat# 13769-1-AP,1:3,000); MMP9 (Proteintech, Cat#10375-2-AP,1:3,000); GAPDH (Proteintech, Cat# 66031-1-1g1:3,000).

Statistical Analysis

SPSS23.0 software (IBM) was performed for statistical analyses. P-values<0.05 were identified as significant. Significant differences were supposed at P<0.05; P<0.01; P<0.001; and P<0.0001.

Results

Identification of Differentially-Expressed IRGs and TFs

We downloaded 351 samples, including 319 OSCC cases and 32 normal cases, and two OSCC samples were excluded owing to insufficient clinical data. In total, 317 OSCC cases and 32 normal controls were enrolled in this study. A total of 24,98 IRGs (Figure 1A and B) were obtained from the ImmPort database and 317 TFs (Figure 1C and D) from Cistrome. Differentially-expressed analysis indicated that there were 381 differentially-expressed IRGs (Supplementary Table S1) and 62 TFs (Supplementary Table S2) in OSCC. According to gene expression level with the cut-off criteria correlation coefficient of more than 0.4, the co-expression network was established in R in order to understand the regulatory relationship between TFs and IRGs. Finally, a total of 25 TFs, 21 IRGs, and 63 edges were enrolled in the co-expression network (Supplementary Table S3), which then was visualized in Cytoscape software (Supplementary Figure S1).

Figure 1 Differentially-expressed analysis of IRGs and TFs. (A, C) The heatmap of IRGs and TFs expression level in TCGA OSCC samples. N represents normal control cases, while T represents tumor cases. The ascending normalized expression level in the heatmaps is colored from blue to red. (B, D) The volcano plot of differentially-expressed IRGs and TFs in OSCC. Red stands for up-regulation and blue stands for down-regulation. Abbreviations: IRGs, immune-related genes; TFs, transcription factors; TCGA, The Cancer Genome Atlas; OSCC, oral squamous cell carcinoma.

Establishment of Gene Expression Signature

Three hundred and eighty-one differentially-expressed IRGs were enrolled in univariate Cox regression analysis, and 42 IRGs were identified as prognostic factors (Figure 2A). Subsequently, 42 prognostic molecules were selected into multivariate Cox regression and step-wise Cox regression analysis. Finally, a gene expression signature including HTN3, CTSG, ADIPOQ, IL1A, CCL26, BTC, IL17F, JAG2, ESRRG, and TNFRSF4 was constructed (Table 1). Risk score=(0.100537503* HTN3)+(−0.155940982* CTSG)+(0.109018512* ADIPOQ)+(0.096145672* IL1A)+(0.099317877* CCL26)+(0.11718248* BTC)+(−0.159415533* IL17F)+(0.167445928* JAG2)+(0.078728308* ESRRG)+(−0.191134673* TNFRSF4). According to the Cox formula median value, each OSCC patient’s risk score was determined and visualized in R software along with corresponding survival time and gene expression profiles (Supplementary Figure S2 A, B). And the low-risk group experienced an excellent overall survival rate compared to the high-risk group (Figure 2B). In addition, CTSG, IL1A, CCL26, and JAG2 were relevant to overall survival (Figure 2C–F), indicating that these four IRGs might be regarded as potential prognostic biomarkers in OSCC.

Table 1 Gene Signature on the Basis of Prognostic IRGs in OSCC

Figure 2 Gene expression signature according to IRGs. (A); Screening prognostic factors through univariate Cox regression analysis. (B); Survival analysis on the basis of immune-related gene expression signature. (C–F) Survival analysis of IRGs enrolled in risk formula.

Stratification Survival Analysis

The prognostic value of each clinical parameter in TCGA was further explored with survival analysis. Distant metastasis was excluded owing to plenty of significantly inadequate data. Survival analyses combined with stratification analyses, including age (Figure 3A), gender (Figure 3B), tumor grade (Figure 3C), clinical stage (Figure 3D), and N (Figure 3E), T (Figure 3F) classification in the TNM system indicated that OSCC patients with a high-risk score had a remarkably undesirable overall survival rate compared to those with low-risk score.

Figure 3 Stratification analysis and survival analysis. Based on gene expression signature, stratification analyses including age (A), gender (B), tumor grade (C), clinical stage (D), N classification (E), and T classification (F) in the TNM system combined with survival analyses were performed in R software “survival” package.

Correlation Analysis

Bioinformatics analyses showed that risk score based on IRGs might be an independent prognostic biomarker in OSCC (Figure 4A and B). Furthermore, correlation analyses indicated that prognostic factors ESRRG and IL1A were also related to T classification in the TNM system (Figure 4C and D). To our surprise, the low CTSG expression level was not only associated with unfavorable overall survival, but also relevant to terminal clinical stage and high T classification (Figure 4E and F). In addition, high-risk score meant higher T classification (Figure 4G) and advanced clinical stage (Figure 4H), and was related to lower levels of immune cells including B-cell (Figure 5A), CD4_T-cell (Figure 5B), CD8_T-cell (Figure 5C), Neutrophil (Figure 5D), Macrophage (Figure 5E), and Dendritic (Figure 5F) in OSCC.

Figure 4 Correlation analysis between risk score and clinical parameters. (A, B) Bioinformatics indicated that risk score may be an independent prognostic biomarker. (C, D) High-risk score meant high T classification in the TNM system and terminal clinical stage (stage Ⅰ+Ⅱ, stage III+Ⅳ). (E, F) High ESRRG and low IL1A expression level meant low T classification. (G, H) Low CTSG expression level meant high T classification in TNM system and terminal clinical stage.

Figure 5 Correlation analysis between risk score and immune cells including B-cell (A), CD4_T-cell (B), CD8_T-cell (C), Neutrophil (D), Macrophage (E), and Dendritic (F) in TCGA OSCC patients.

CTSG Inhibits OSCC Cells Proliferation, Migration/Invasion

CTSG, a member of the peptidase S1 protein family, is found in azurophil granules of neutrophilic polymorphonuclear leukocytes. It may participate in the killing and digestion of engulfed pathogens, and connective tissue remodeling at sites of inflammation. In our study, bioinformatics analyses indicated that CTSG may be an independent prognostic biomarker in OSCC. Therefore, CTSG expression level was detected in OSCC (Figure 6A and B) tissues. Western blot assay showed that CTSG protein level was downregulated in eight out of 10 OSCC cases (Figure 6C). Subsequently, CTSG overexpression vectors were transfected into SCC9 and SCC25 cell lines, and the transfection efficiency was detected by QPCR and Western blot (Figure 7A and B). As shown in the results, overexpression of CTSG inhibited the EMT pathway, which is remarkably relevant to cancer initiation and progression¹⁸ in OSCC via upregulation of E-cadherin and downregulation of N-cadherin, MMP2, MMP9, and Vimentin in SCC9 and SCC25 cells (Figure 7C). The protein–protein interaction network also indicated that CTSG and MMP9 had interactive regulation (SupplementaryFigure S2C). Finally, the results of CCK8, clonogenic, and transwell assay indicated that overexpression of CTSG obviously inhibited the proliferation, migration, and invasion in OSCC cells (Figure 7D and E; Figure 8A–E).

Figure 6 CTSG was differentially-expressed in OSCC. (A) CTSG was lower expression level in SCC9 (P=0.0021); SCC15 (P=0.0018); SCC25 (P=0.0012). (B) CTSG relative CT value in OSCC specimens and adjacent normal specimens (P=0.0433). The higher the relative CT value, the lower the CTSG expression level. (C) CTSG protein level in 10 OSCC cases. Eight cases were lower protein level in OSCC tumor tissues. N meant adjacent normal control, while T stands for tumor tissues. Significant differences were supposed at *P<0.05; **P<0.01.

Figure 7 CTSG inhibited the activity of EMT pathway and OSCC progression. (A, B) Transfection efficiency was detected by qRT-PCR and Western blot assay (SCC9, P<0.0001; SCC25, P=0.0003). (C) Overexpression of CTSG up-regulated E-cadherin and down-regulated N-cadherin, MMP2, MMP9, and Vimentin both in SCC9 and SCC25 cell lines. (D) SCC9 cell lines clonogenicity (P=0.0103), migration (P=0.0017), and invasion (P=0.0012) statistical analyses. (E), SCC25 cell lines clonogenicity (P=0.0236), migration (P=0.0003), and invasion (P=0.0021) statistical analyses. NC stands for negative control, while CTSG-OE represents CTSG overexpression. Significant differences were supposed at *P<0.05; **P<0.01; ***P<0.001; ****P<0.0001.

Figure 8 CTSG inhibits OSCC cells proliferation, migration, and invasion. (A, B) CCK8 assay indicated that high CTSG expression inhibits OSCC cells proliferation (SCC9; 3-day, P=0.0018, 4 days, P=0.0045; SCC25; 3-day, P=0.0235, 4 days, P=0.0067). (C) CTSG overexpression inhibited the clonogenicity in vitro. (D, E) CTSG overexpression reduced the migration and invasion ability of OSCC cell lines. NC stands for negative control, while CTSG-OE represents CTSG overexpression. Significant differences were supposed at *P<0.05; **P<0.01.

Discussion

Owing to the high recurrence rate and unsatisfied prognosis, OSCC poses a great challenge to modern healthcare. Recently, studies have demonstrated that the immune system was associated with OSCC initiation and progression.⁷ Therefore, it might be helpful to identify effective biomarkers of OSCC through exploring the role of whole IRGs by integrated bioinformatics analysis.

In our study, 381 IRGs and 62 TFs were identified to be differentially expressed in OSCC. A co-expression network was constructed successfully to investigate the regulatory relationship between IRGs and TFs. Some node genes were associated with OSCC. For instance, FOXM1 was identified as a novel predictor of tumor recurrence and its potential involvement in epithelial–mesenchymal transition (EMT) in OSCC.¹⁹ BIRC5, regulated by FOXM1 in the co-expression network, was an independent prognostic factor for patients with OSCC.²⁰ In addition, E2F1, a transcription factor of BIRC5, was also relevant to OSCC progression.²¹ HIF1A regulated RGMA expression which was associated with OSCC proliferation.²² TGFB1, regulated by HIF1A in the co-expression network, could induce THBS1 to promote OSCC invasion.²³ These node genes played a vital role in OSCC initiation and progression; therefore, the specific regulatory relationship between these TFs and IRGs should be further explored. Subsequently, a gene expression signature according to prognostic IRGs was built which could predict the prognosis of OSCC accurately. Correlation analyses indicated that risk level was related to T classification and clinical stage. In addition, the higher risk score meant the lower immune cells level and terrible overall survival in OSCC patients, which indicated that the immune function deficiency was closely relevant to the prognosis of patients with OSCC. Integrated bioinformatics analysis indicated that risk level according to IRGs might be an independent prognostic biomarker. Therefore, the 10 IRGs may be prognostic biomarkers of OSCC. Most of them were proved to be related to various cancers. For example, CTSG could promote cell migration and multicellular aggregation in human breast cancer MCF-7 cells through depending on E-cadherin²⁴ and was an immunotherapeutic target in acute myeloid leukemia (AML).²⁵ Our study showed that CTSG could inhibit OSCC proliferation, migration, invasion, and the induction of EMT in vitro, which plays a key role in tumorigenesis and progression, and a higher CTSG expression level meant well-pleasing overall survival and early-stage tumor. Our results indicated that CTSG plays a vital inhibitory role in OSCC carcinogenesis and may be a potential biomarker and therapeutic target for OSCC. The study demonstrated that CTSG could activate MMP2 which was associated with inflammatory cells in tumor invasion and angiogenesis. According to the CTSG protein–protein interaction network, CTSG, MMP9, and SERPINB4 had a mutual regulation relationship. Recently, numerous studies indicated that MMP9 might be involved in tumor invasion, metastasis,²⁶ angiogenesis, and local tumor progression in OSCC,²⁷ which indicating that CTSG might play a significant role in OSCC by interaction with MMP9. In addition, Masashi Shiiba et al²⁸ suggested that SERPIN-family genes, including SERPINB4, played a crucial role in OSCC carcinogenesis and might be involved in tumor suppression, which showed that CTSG might have an interaction with these genes and a close relation with OSCC initiation and progression. The study showed that extracellular mast cell-derived IL-17F at the tumor invasion front was associated with better disease-specific survival in patients with all-stages and early-stages of oral tongue SCC²⁹ and serum IL-17F combined with vascular endothelial growth factor (VEGF) served as potential diagnostic biomarkers for OSCC.³⁰ Furthermore, Hatano et al³¹ identified that overexpression of JAG2 was associated with undesirable overall survival, and JAG2 may be a potential prognostic biomarker in OSCC, which was similar to our results. ESRRG was verified as a negative regulator of the Wnt signaling pathway in gastric cancer (GC).³² Li et al³³ demonstrated that IL1A was related to the proliferation of cervical carcinoma cells. In our study, these IRGs were closely correlated with the prognosis of OSCC and might be potential biomarkers in OSCC via bioinformatic analysis. However, the roles of them in OSCC remain obscure, and further studies should be made to explore their mechanism.

Conclusions

In conclusion, our studies showed that IRGs were involved in OSCC initiation, progression, and prognosis via bioinformatics analyses. A co-expression network was constructed to predict the regulatory relationship between IRGs and TFs. In addition, a gene expression signature according to prognostic IRGs could predict OSCC patients' prognosis and immune cells level accurately, and the risk-score was served as an independent prognostic biomarker in OSCC, indicating that these IRGs might play a crucial role in OSCC progression. Finally, CTSG played an important role in OSCC carcinogenesis and might be a potential immunotherapeutic target for OSCC.