CBX3 is a Prognostic Biomarker Correlated with ATR Activation and Immune Infiltration in Head and Neck Squamous Cell Carcinoma

Introduction

Head and neck cancer is the sixth leading malignancy worldwide, and nearly 90% of cases are classified as head and neck squamous cell carcinoma (HNSCC).¹ HNSCC is a heterogeneous solid tumor with an aggressive phenotype and poor clinical outcomes, primarily due to local recurrence and regional lymph-node and distant metastasis. With the vast improvement in treatments, including surgery, radiotherapy, chemotherapy, and multidisciplinary comprehensive sequence therapy, the quality of life of patients with HNSCC has improved to a certain extent; however, the five-year survival rate remains at approximately 50%.² Although the knowledge of tumor invasion and metastasis has improved, the molecular mechanism of HNSCC tumorigenesis remains unclear.

In humans, eight chromobox (CBX) proteins have been identified, with links primarily to heterochromatin, gene expression, apoptosis, and the regulation of development.³ CBX family members are abnormally overexpressed in tumorigenesis and affect tumor prognosis.^4,5 For example, overexpression of CBX8 in hepatocellular carcinoma patients is related to short overall survival (OS) time.⁶ In gastric cancer, CBX7 is activated via the AKT pathway and downregulated by the tumor suppressor protein p16, and plays an important role in maintaining the stem cell-like characteristics of gastric cancer cells.⁷ However, the roles of the various CBX family members in the onset and development of HNSCC remain unclear.

In the present study, we used a range of tools and databases to investigate the relationships between CBX family members and HNSCC. We analyzed the expression, mutation, and clinical correlations of CBX proteins in HNSCC patients to predict their potential functions and prognostic values.

Materials and Methods

Data Collection and Bioinformatics Analysis

To analyze the expression and correlations with HNSCC prognosis of CBX family members, we downloaded gene expression data (RNA-sequencing [RNA-seq] data expressed as transcripts per million) from The Cancer Genome Atlas (TCGA) (http://cancergenome.nih.gov). In total, 546 HNSCC cases were included in the analysis, for which RNA-seq data and clinicopathological data were extracted. Patient characteristics included sex, race, TNM stage, clinical stage, histologic grade, lymphovascular invasion, smoking status, and radiation therapy or not. All data used in this study were obtained from TCGA, thus, ethics approval and informed consent were not required, which was approved by the Ethics Committee of the Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine. We also collected 17 pairs of HNSCC tumor and normal tissues for.

RNA-Seq Gene Expression Analysis

The EdgeR package (version 3.8, http://www.bioconductor.org/EdgeR) was used to analyze the differentially expressed mRNAs between HNSCC tissues and normal adjacent tissues.⁸ A |log(fold change)| > 2 and P < 0.01 were set as cut-off criteria. RNA expression was visualized in a histogram using the ggplot2 package (version 3.1.0, https://github.com/tidyverse/ggplot2).

Survival Prognosis Analysis

We used the survminer package to obtain OS and disease-free survival significance map data of the CBX family across all TCGA tumors. Cutoff-high (50%) and cutoff-low (50%) values were used as expression thresholds for splitting the high-expression and low-expression cohorts. The Log rank test was used in the hypothesis test, and survival plots were constructed using the survminer package. The relationships between clinicopathologic features and CBX expression were analyzed using the Wilcoxon rank-sum test and logistic regression. Clinicopathologic characteristics associated with OS and DSS were analyzed with Cox regression and the Kaplan-Meier method. The multivariate Cox analysis was used to identify the influence of CBX expression on survival, along with other clinical features. All hypothetical tests were two-sided, and P values<0.05 were considered significant.

Functional Enrichment Analysis

Metascape (http://metascape.org) is a free gene-list analysis tool for gene annotation and analysis.⁹ We used Metascape to conduct pathway and process enrichment analysis of differentially expressed genes (DEGS) in HNSCC significantly associated with CBX3 alterations. Gene enrichment was analyzed using Gene Ontology (GO) terms for biological process, cellular component, and molecular function categories, as well as Kyoto Encyclopedia of Gene and Genomes (KEGG) pathways. Terms with P < 0.01, minimum count of 3, and enrichment factor > 1.5 were considered as significant.

Gene Set Enrichment Analysis (GSEA)

We used GSEA to list all DEGs according to their correlation with CBX3 expression. The R package clusterProfiler (version 3.6.3) was used to perform GSEA of the high- and low-CBX3 groups. For each analysis, gene set permutation was repeated 1000 times according to the default statistical method. DEGs with a false discovery rate-adjusted P-value (q-value) < 0.25, adjusted P < 0.05, and |normalized enrichment score| > 1 were considered as significantly enriched.

Analysis of Immune Cell Characteristics by Single-Sample (ss)GSEA

We used the GSVA R package to compare the normalized CBX3 gene expression data and immune cell signatures. ssGSEA classifies marker gene sets with common biological functions, chromosomal localization, and physiological regulation in a single sample. Immune infiltration in HNSCC was analyzed using ssGSEA in the R package (version 3.6.3). To reveal correlations between CBX3 and immune cell infiltration levels, we used the Spearman’s rank correlation test and Wilcoxon’s rank-sum test (P < 0.05).

Tissue Samples and qRT-PCR

Paired, fresh HNSCC tissues and their corresponding adjacent normal tissues were collected from the Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine. Written informed consent was obtained from each patient. The patients were well-informed and have consented to this study, and the process was approved by Ethics Committee of the Ninth People’s Hospital, Shanghai Jiao Tong University School of Medicine.(SH9H-2021-TK547-1).

Total RNA was extracted using TRIzol (Invitrogen, CA, United States). An equal amount of RNA was reverse-transcribed using the HiScript II Q RT Supermix and was quantified by qPCR using SYBR Green (Bimake). The primer sequences are shown in Supplementary Table 1.

Results

CBX1–8 Expression in Different Cancers

The public database, TCGA, was used to analyze the mRNA expression of CBX1–8 in various cancers and normal tissues. Figure 1 shows the expression levels in pan-cancer. CBX1, CBX2, CBX3, CBX4, CBX5, and CBX8 mRNA levels were elevated in most cancers, whereas CBX6 and CBX7 were downregulated in most cancers.

Figure 1 Expression level of CBX family was upregulated in most cancer types.

CBX1–8 Expression Differs Between HNSCC Tissues and Normal Tissues

Next, we focused on the expression of CBX1–8 in HNSCC using TCGA data. CBX1, CBX2, CBX3, CBX4, CBX5, CBX6, and CBX8 mRNA levels were higher in tumor tissues than in normal tissues, whereas that of CBX7 was lower in tumor tissues (Figure 2A). Receiver operating characteristic (ROC) curve analysis revealed that most CBX family members, including CBX1, CBX2, CBX3, CBX8, may be used as a marker to distinguish tumor from non-tumor cells (Figure 2B). The cut-off levels of CBX1-8 are shown in Supplementary Table 2. We also used qPCR to validate mRNA expression levels of CBX3. The results showed that CBX3 mRNA levels were higher in tumor tissues than in normal tissues (Figure 2C).

Figure 2 Expression level of most members of CBX family was upregulated in HNSCC. (A) CBX expression of normal samples and tumor samples in TCGA head and neck squamous cell carcinoma. (B) ROC curve of CBX family. (C) mRNA relative expression of CBX3 (qPCR results).

Prognostic Value of CBX1–8 Expression in HNSCC

Multivariate Cox regression of OS and disease-specific survival (DSS) showed that CBX3 and CBX5 expression levels were independent risk factors in patients with HNSCC (Table 1A and B). Therefore, we focused on these two genes and used Kaplan–Meier curves of OS and DSS to analyze the prognosis of HNSCC patients with different expression levels of CBX3 and CBX5. High CBX3 expression was significantly associated with shorter survival, while high CBX5 expression was associated with longer survival (Figure 3A and B). We used a subgroup analysis to analyze the relationships between CBX3 and CBX5 expression and the prognosis of patients with different characteristics (Figure 3C and D). In patients undergoing radiation-therapy, patients with clinical grades III and IV, patients with tumor grades III and IV, and smoking patients, high CBX3 expression predicted poor survival, whereas high CBX5 expression predicted good survival in these patient groups.

Table 1 (A) CBX3 Expression Associated with Clinicopathological Parameters (Logistic Regression). (B) CBX5 Expression Associated with Clinicopathological Parameters (Logistic Regression)

Figure 3 Prognostic value of CBX3 and CBX5 in HNSCC. (A) Kaplan-Meier curves of OS and DSS with different expression level of CBX3. (B) Kaplan-Meier curves of OS and DSS with different expression level of CBX5. (C) Kaplan-Meier curves of OS and DSS with different expression level of CBX3 in radiation therapy group, clinical grades III and IV group, T grades III & IV group and smoked group. (D) Kaplan-Meier curves of OS and DSS with different expression level of CBX5 in radiation therapy group, clinical grades III and IV group, T grades III & IV group and smoked group.

CBX3 and CBX5 Expression Correlate with Clinicopathological Parameters

We analyzed the correlation of CBX3 and CBX5 expression with clinicopathological parameters. We used clinical information of 502 HNSCC patients, including 368 male and 124 female patients with a median age of 61 years (range: 55–71 years), collected from the TCGA. Clinicopathological features are shown in Table 2. CBX3 and CBX5 expression in HNSCC tissues was categorized as “low” or “high” based on median values. High CBX3 expression was strongly associated with high T&N stages, age, clinical stages, radiation therapy, high histologic grade, and lymphovascular invasion (Table 3A and B), while high CBX5 expression was strongly associated with high histologic grade.

Table 2 Characteristic of Patients with HNSCC Based on TCGA

Table 3 (A) Univariate and Multivariate Cox Regression Analysis of Clinicopathological Parameters and Overall Survival. (B) Univariate and Multivariate Cox Regression Analysis of Clinicopathological Parameters and Disease-Specific Survival

Identification of DEGs in High- and Low-CBX3 Expression Samples via Functional Cluster Analysis

To explore the potential mechanism of CBX3 in causing tumor progression, we analyzed DEGs in high- and low-CBX3 expression samples. In total, we identified 273 DEGs, including 85 upregulated and 188 downregulated genes. DEG expression is presented in a heatmap and a volcano plot (Figure 4A and B). The functions of the DEGs were analyzed by GO and KEGG enrichment analyses. The top GO terms in the biological process, molecular function, and cellular component categories were cell adhesion molecule binding, cornified envelope, and peptidase regulator activity, respectively. KEGG analysis showed that genes involved in Staphylococcus aureus infection were enriched (Figure 4C). We used GSEA to analyze the key pathways related to CBX3. The most differentially enriched pathways based on the normalized enrichment scores were 1) activation of ataxia-telangiectasia mutated and Rad3-related protein kinase (ATR) in response to replication stress and 2) activation of the pre-replicative complex (Figure 4D).

Figure 4 Differentially expressed genes (DEGs) in High CBX3 expression and Low CBX3 expression samples and functional cluster analysis. (A) Volcano map of DEGs with |log2FoldChange >2| and adjusted P value<0.05. (B) Heat map of top 10 DEGs in high and low CBX3 expression groups. (C) Significantly enriched GO and KEGG annotations of CBX3 related genes. (D) Enrichment genes analysis by GSEA.

Correlation Between CBX3 Expression and Immune Infiltration

We used Spearman correlation to analyze the correlation between CBX3 expression and immune cell infiltration levels (generated by ssGSEA). The results are shown in Figure 5A. CBX3 expression was significantly negatively correlated with mast cells. The Spearman R value was –0.314, with P < 0.001. The Wilcoxon rank-sum test corroborated that compared with the CBX3-high expression group, mast cell infiltration was higher in the CBX3-low expression group (Figure 5B). Other immune cells that were statistically significantly enriched were dendritic cells (DCs), neutrophils, immature DCs, and Th17 cells (Figure 5B).

Figure 5 The expression level of CBX3 was associated with the immune infiltration in tumor environment. (A) The plots showed the correlation between CBX3 expression and immune cells subsets. (B) Spearman correlation and expression distribution of mast cells, DCs, neutrophils, iDCs, and Th17 cells in low and high CBX3 samples.

Discussion

Patients with HNSCC have a poor five-year survival (<50%) and a high risk for regional recurrence or distal metastasis.¹⁰ Efficient biomarkers for HNSCC diagnosis and treatment are lacking. The identification of novel biomarkers for HNSCC diagnosis and prediction is needed. CBX3 interacts with various cofactors to exert various biological functions, including RNA alternative splicing, DNA damage response, transcription elongation, cell growth, and differentiation.¹¹ CBX3 is normally not detected in differentiated cells of various normal human tissues, but is detected in tumor tissues of various cancers, including esophageal, lung, cervical, colon, and breast cancers.^12–15 Substantial evidence indicates that CBX3 has cancerogenic and prognostic roles in various cancer types. However, the exact functions of CBX3 in HNSCC remain unclear.

In this study, by analyzing TCGA data, we found that CBX1, CBX2, CBX3, CBX4, CBX5, and CBX8 exhibit increased expression in tumor tissues of most cancers when compared with normal tissues, while CBX6 and CBX7 show reduced expression. In HNSCC, many CBX family members were highly expressed in tumor tissues. The ROC curve showed that CBX1, CBX2, CBX3, CBX4, CBX5, and CBX8 exhibited potential as diagnostic makers.

We researched the relationships between CBX family members and prognosis. High CBX3 expression was correlated with poor prognosis, whereas high CBX5 expression was correlated with good prognosis. The other members were not correlated with prognosis. Therefore, we further focused on CBX3 and CBX5. The odds ratios of CBX3 and CBX5 expression indicated that they may serve as diagnostic markers. Moreover, high CBX3 expression in HNSCC was associated with a high clinical stage, high histologic grade, and lymphovascular invasion. High CBX5 expression in HNSCC was associated with high histologic grade. CBX3 performed well as a prognostic marker as indicated by these results.

To investigate the potential function of CBX3 in HNSCC further, we analyzed DEGs in low- and high-CBX3 expression groups. Among the upregulated DEGs, AIMP2, SRSF2, and ERH are reportedly related to tumor progression in several cancers. High AIMP2 expression predicts poor prognosis in hematologic cancer.¹⁶ SRSF2 is a key regulator of RNA splicing dysregulation in cancer with possible clinical implications as a candidate prognostic factor in patients with hepatocellular carcinoma.¹⁷ ERH regulates MYC gene expression and is related to the migration and invasion of bladder cancer.¹⁸ According to GO analysis, CBX3 affects adhesion molecule binding and epidermis development, which is associated with epithelial-mesenchymal transition of cancer cells.¹⁹ According to KEGG analysis, CBX3 involved in Staphylococcus aureus infection pathway, but there was no research about the relationship between CBX3and this pathway. These results suggest that CBX3 may affect cancer progression via various mechanisms.

We performed GSEA to search for potential pathways affected by CBX3. The results indicated that CBX3 may affect HNSCC progression through ATR activation in response to replication stress and activation of the pre-replicative complex, which are related to DNA repair in cancer cells.²⁰ ATR is a critical component of the cellular DNA damage response.²¹ ATR is activated by regions of single-stranded DNA, some of which occur as a result of replication stress.^22–24 High ATR expression may provide a stable environment for oncogenes and relate to poor prognosis and drug resistance.²⁵

Tumor-infiltrating immune cells offer prognostic value in several human malignancies.^26–28 Therefore, we analyzed the correlations between CBX3 expression and immune cells. CBX3 expression was negatively associated with mast cells, DCs, immature DCs, and neutrophils. These immune cells can act as antigen-presenting cells (APCs) that ingest, process, and present extracellular antigens, activate CD4+ T cells, and induce an immune response.²⁸ High CBX3 expression was associated with a reduction APCs in this study, which implies that CBX3 may suppress adaptive immunity and killer T-cell function by inhibiting the function and proliferation of APCs. However, the detailed underlying mechanism requires further investigation.

Conclusions

In conclusion, our research showed that high CBX3 expression predicts decreased OS and that CBX3 can function as a potential prognostic marker. CBX3 may affect cancer progression via various mechanisms that remain to be elucidated and promote HNSCC proliferation and invasion by activating ATR and immune cell infiltration.