ONCOMINE (www.oncomine.org) is an accessible online cancer microarray database providing powerful, genome-wide expression analysis (Rhodes et al., 2004). In our study, data were obtained to assess the transcriptional expression of CBX proteins between different carcinomas and adjacent normal control tissues ESCA. The difference in transcriptional expression of CBXs in ESCA was determined by Student’s t-test. Sufficient fold changes ≥1.5, significant p value < 0.05, and gene rank ≥ the top 10% were set as the threshold.

We obtained the raw counts of RNA-sequencing data and corresponding clinical information of CBX family members of TCGA ESCA tissue data from the TCGA dataset (http://tcga-data.nic.nih.gov/). RNAseq data in FPKM (Fragments Per kilobase per Million) format and log2 conversion for expression comparison between samples (Vivian et al., 2017). We draw receiver-operating characteristic (ROC) curves by using the “pROC” package. For Kaplan–Meier curves, p values and the hazard ratio (HR) with 95% confidence interval (CI) were generated by log-rank tests and unvaried Cox proportional hazards regression. CBX family member expressions and their correlation with the infiltration abundance of immune cells such as ADC [activated DC], B cells, CD8 T cells, Cytotoxic cells, DC, Eosinophils, and iDC [immature DC] in ESCA were evaluated using Spearman’s correlation with TCGA ESCA in the project level 3 HTSeq-RNAseq FPKM format data and clinical data. All analytical methods previously mentioned and R packages were performed using R software version v3.3.6 (Liu et al., 2018). p < 0.05 was regarded as statistically significant.

The gene expression microarray of GSE20347 and GSE38129 were downloaded from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) of the National Center for Biotechnology Information (NCBI) via the GEO 1010query package (Davis and Meltzer, 2007). Altogether, 17 tumors and matched normal adjacent tissue samples were obtained from GSE20347, while 30 ESCC and 30 normal samples were obtained from GSE38129. The detection platform of the above expression microarrays was GPL571 [HG-U133A_2] Affymetrix Human Genome U133A 2.0 Array.

GEPIA (http://gepia.cancer-pku.cn/) is a newly generated web server containing RNA sequence expression data of 9,736 tumors and 8,587 normal samples based on TCGA and the GTEx databases (Tang et al., 2017). We used the “Single Gene Analysis” module to perform differential mRNA expression analysis according to pathological stages, survival analysis, and correlation analysis. “Multiple Gene Comparison” module of GEPIA was utilized to assess the multiple gene comparison analysis of the CBX family, using the “ESCA” dataset. The p value cutoff was 0.05.

UALCAN (http://ualcan.path.uab.edu/analysis.html) is a comprehensive database based on level 3 RNA-seq and clinical data of 31 cancer types from The Cancer Genome Atlas (TCGA) and MET500 cohort data (Chandrashekar et al., 2017). Therefore, it can provide information about the relative transcriptional expression of genes in carcinomas compared with normal samples. Furthermore, the database also presents information on the association of transcriptional expression with relative clinicopathological features. In our study, we employed UALCAN to investigate the CBXs mRNA expressions and their relationship with clinicopathological parameters of ESCA. The significant difference of transcriptional expression between groups was evaluated using two-sample Student’s t-test and p < 0.01 was regarded as statically significant.

The Kaplan–Meier plotter (http://kmplot.com/analysis/) is a tool that can give information on the effects of 54,000 genes on survival in 21 cancer types (Nagy et al., 2018). We made use of the Kaplan–Meier plotter to evaluate the association between the mRNA expression levels of CBX members with OS of ESCA patients. Information about the number-at-risk cases, median values of mRNA expression levels, the hazard ratio (HR), 95% confidence intervals (CIs), and log-rank p value can be accessed at the K-M plotter webpage. A statically significant difference was considered when p value was below 0.05.

cBioPortal (www.cbioportal.org) is an open-access resource used to visualize and analyze multidimensional cancer genomics datasets (Cerami et al., 2012). Based on the TCGA database, genetic alterations of the CBX gene in patients with ESCA were obtained from cBioPortal. The obtained mRNA expression z-score threshold was 1.8 between the unaltered and altered patients.

Metascape (http://metascape.org) is a free, well-maintained, user-friendly tool for gene annotation and gene list enrichment analysis (Zhou et al., 2019). We used Metascape to perform pathway and process enrichment analysis of the CBX genes and neighboring genes closely related to CBXs alteration via the “Custom Analysis” module, for GO and KEGG enrichment as well as protein-protein interaction analyses.

TIMER (https://cistrome.shinyapps.io/timer) is a comprehensive resource that could provide systematic analyses with the dataset of 10,897 samples among diverse cancers in the TCGA database (Li et al., 2017). CBX family expression scatter plots and their correlation with the abundance of immune cells such as B cells, CD4⁺ T cells, CD8⁺ T cells, neutrophils, macrophages, and dendritic cells in ESCA were assessed using Spearman’s correlation with TCGA_ESCA datasets. The infiltration abundance for each somatic copy number alterations (SCNA) category was compared to the normal by a two-sided Wilcoxon rank-sum test, and statistical significance was identified as p < 0.01.

ESCA tissues and the corresponding adjacent normal tissues were collected from patients undergoing surgery at Fujian Union Hospital between April 2020 and April 2021. Patients who received neoadjuvant chemoradiotherapy were excluded. The Ethics Committee of the Fujian Medical University Union Hospital approved for use of all specimens, and all patients provided written informed consent. The tissues were stored in liquid nitrogen until use.

We quantified the level of CBX family members in tumor tissues (n = 17) as well as adjacent normal tissues (n = 17) obtained from patients with ESCA, using Quantitative real-time PCR. Total RNA from specimens were isolated using Trizol reagent (Invitrogen). Quantitative RT-PCR was performed with SYBR Green Real-Time Mix (Roche) by a 7,500 Real-time PCR according to the manufacturer’s protocol (Applied Biosystems). PCR amplification was conducted in the following conditions: 95°C for 10 min, 95°C for 15 s (denature), and 60 °C for 1 min (anneal/extend) for 40 cycles, 95 °C for 15 s, and 60 °C for 1 min, and then 95 °C for 15 s (Mel curve). Primers are listed in Supplementary Table S1. The fold change was quantified via 2 − ΔΔCt [ΔΔCt = (ΔCt of genes of interest) − (ΔCt of β-actin)].

We extracted protein from tumor tissues (n = 12) and adjacent normal tissues (n = 12) tissues in radioimmunoprecipitation assay (RIPA) buffer (Beyotime, Shanghai, China). The protein content was determined by the Bicinchoninic acid (BCA) assay kit (Beyotime, Shanghai, China). Proteins were separated by SDS-PAGE and were visualized using the ECL system (Bio-Rad, Hercules, CA, United States). Primary antibodies used in this study were as follows: CBX1 Polyclonal Antibody (Proteintech, 10241-2 -AP), CBX3 Polyclonal Antibody (Proteintech, 11650-2 -AP), CBX7 Polyclonal Antibody (Proteintech, 26278-1-AP), and CBX8 Monoclonal Antibody (Santa Cruz, sc-374332). Monoclonal Anti-GAPDH antibody was purchased from Cell Signaling Technology (Danvers, MA, United States).

ESCA tumor tissues and adjacent normal tissue were collected after surgery, then were fixed in 10% neutral buffered formalin and embedded in paraffin. Tissue samples were then sectioned at 4um thickness. Sections were deparaffinized, rehydrated, processed for antigen retrieval, blocked, incubated with primary antibody (CBX1:1:500 dilution, CBX3:1:500 dilution, CBX8:1:50 dilution) at 4 °C overnight, followed by incubation corresponding fluorescence-conjugated secondary antibody (1:500 dilution) for 1 h at room temperature. Then slides were incubated with DAPI (BD Biosciences) for 10 min at RT and mounted with ProLong Gold antifade reagent (Invitrogen). Subsequently, the fluorescence images were captured with confocal microscopy.

The expression of CBX family members were collected using ONCOMINE database (Figure 1A). We first investigated the mRNA transcriptional levels of CBXs in ESCA and normal tissues with ONCOMINE. As shown in Figure 1A and Table 1, the transcriptional levels of CBX1, CBX3, and CBX6 in ESCA were significantly elevated while CBX7 was significantly reduced compared to normal tissues. TCGA ESCA cohort showed that the expression levels of CBX1 (p < 0.001), CBX2 (p < 0.001), CBX3 (p < 0.001), CBX4 (p < 0.001), and CBX8 (p < 0.001) were higher while the expression level of CBX7 (p < 0.001) was lower in ESCA samples compared to paired para-cancerous samples (Figure 1B). Then we evaluated the expression levels of CBXs in ESCA and normal tissues with the GEO database (GSE20347 and GSE38129). The transcriptional levels of CBX1 (p < 0.001), CBX3 (p < 0.001), CBX5 (p < 0.001), CBX6 (p = 0.006), and CBX8 (p = 0.049) in ESCA were significantly elevated, while CBX7 (p = 0.015) was decreased in ESCA tissues compared to normal tissues (Figure 1C). We also compared the relative expression levels of CBXs in ESCA using GEPIA, the results showed that the relative expression of CBX3 was the highest among all the CBX proteins (Figure 1D).

Next, we used the ROC curve of CBX family members to access the diagnostic capability of CBX family members for ESCA by using the TCGA cohort. The results indicated that CBX1, CBX2, CBX3, CBX4 and CBX8 had great diagnostic capability with AUC of 0.916(95% CI: 0.838–0.994),0.911(95% CI: 0.837–0.986),0.926(95% CI: 0.841–1.000),0.946(95% CI: 0.872–1.000), and 0.911 (95% CI: 0.820–1.000), respectively. CBX7 had moderate diagnostic capability with AUC of 0.788 (95% CI: 0.681–0.896). CBX5 and CBX6 had low diagnostic capability with AUC of 0.616 (95% CI: 0.451–0.780) and 0.515 (95% CI: 0.381–0.790) (Figure 2), respectively.

In short, we found CBX1, CBX3, and CBX8 were significantly upregulated in ESCA vs. normal tissues, while CBX7 was downregulated in ESCA. Furthermore, CBX1, CBX2, CBX3, CBX4, and CBX8 had great diagnostic capability to distinguish ESCA from normal tissues.

Since the mRNA expression of CBX family members was aberrant in ESCA patients, we next analyzed the relationship between the mRNA expression of different CBX family members and the clinicopathological parameters of ESCA patients. First, we utilized GEPIA to further detect the correlation between the mRNA expression levels of different CBX family members and the pathological stage of ESCA patients. We found a significant correlation between the expression of CBX1 (p = 0.004) and the pathological stage (Figure 3A). CBX1 expression was increased in the advanced stage (stage II/III) as compared with those in the early tumor stage (stage I). These data suggest that CBX1 plays a significant role in the tumorigenesis and progression of ESCA.

We also analyzed the relationship between the mRNA expression of different CBX family members and the cancer stages, tumor grades, tumor nodal metastasis status, and TP53 mutation status of ESCA patients by UALCAN (http://ualcan.path.uab.edu). As shown in Figure 3B, mRNA expressions of CBX1, CBX2, CBX3, CBX4, and CBX8 were upregulated in ESCA as compared with normal. The expression of CBX1 tended to be higher as the tumor stage increased, which was consistent with the previous findings in the GEPIA database. Nevertheless, the mRNA expression of CBX7 was the highest in normal tissues, and as the cancer stage increased, the mRNA expression of CBX7 tended to decrease. CBX6 mRNA expression had no significant relationship with tumor stages. Similarly, mRNA expressions of CBX1, CBX2, CBX3, CBX4, and CBX8 were significantly related to patients’ tumor nodal metastasis status. Patients who were in more advanced tumor nodal metastasis status tended to express higher mRNA of CBX1. However, normal tissues had the highest mRNA expression of CBX7, and the mRNA expression of CBX7 tended to be lower as tumor nodal metastasis status increased (Figure 3C). We then assessed the relationship between the different expressions of CBX family members and the tumor grades of ESCA patients. As shown in Supplementary Figure S1A, mRNA expressions of CBX1, CBX2, CBX3, CBX4, and CBX8 were significantly associated with tumor grades. The mRNA expression of CBX1, CBX2, CBX3, CBX4, and CBX8 tended to be elevated as the tumor grade increased. However, the mRNA expression of CBX7 tended to be lower in grade 4 compared with grade 3. TP53 mutation has been recognized as the most common event and frequently occurs in ESCA (Yang et al., 2020). Interestingly, as shown in Supplementary Figure S1B, CBX family members’ expressions were associated with TP53 mutation status in ESCA patients. CBX1, CBX2, CBX3, and CBX8 were upregulated in ESCA patients with TP53 mutation compared with normal tissues and TP53 non-mutation patients.

Taken together, these data suggested that the mRNA overexpression of CBX1, CBX3, and CBX8 were significantly related to tumor grades and patients’ tumor nodal metastasis status. CBX1, CBX2, CBX3, and CBX8 were significantly related to TP53 mutation in ESCA patients.

To evaluate the value of differentially expressed CBX family members in the progression of ESCA, we assessed the correlation between differentially expressed CBX family members and clinical outcomes using GEPIA. Disease-free survival curves (DFS) were presented in Supplementary Figure S2A. ESCA patients with lower transcriptional levels of CBX1 (HR = 1.6, p = 0.044) were significantly associated with longer DFS. The value of differentially expressed CBX in the overall survival of ESCA patients was also evaluated (Supplementary Figure S2B). CBX family members did not seem to have a significant effect on Overall Survival (OS).

We also analyzed the prognostic values of CBX family members in patients with ESCA using the TCGA database. Also, the result suggested that ESCA patients with higher CBX3 ((HR = 1.78, p = 0.026)) and CBX4 (HR = 1.94, p = 0.011) levels had shorter OS (Supplementary Figure S3A). Increased CBX4 (HR = 2.07, p = 0.017) and CBX5 (HR = 1.94, p = 0.03) mRNA levels were remarkably associated with shorter DFS (Supplementary Figure S3B). Moreover, increased CBX 1 (HR = 1.71, p = 0.045) mRNA level was significantly correlated with shorter Progress Free Interval (PFI) (Supplementary Figure S3C).

Likewise, we used the Kaplan–Meier plotter to analyze the prognostic values of CBX family members in patients with ESCA subgroup analysis (Figure 4). Esophageal squamous cell carcinoma patients with the higher transcriptional level of CBX4 (HR = 2.93, p = 0.008) were closely associated with shorter OS (Figure Figure4A). The overexpression of CBX3 (HR = 3.12, p = 0.00028) and CBX8 (HR = 2.27, p = 0.035) mRNAs in esophageal adenocarcinoma patients were significantly correlated with shorter OS, whereas the overexpression of CBX7 mRNA (HR = 0.48, p = 0.039) was significantly correlated with longer OS (Figure 4B).

Overall, these data demonstrated that increased CBX1/3/4/5/8 and decreased CBX7 were significantly associated with unfavorable clinical outcomes in ESCA patients.

Using cBioPortal, we investigated the genetic alterations of CBX family members and discovered a high alteration frequency (52%) in ESCA patients. CBX3, CBX1, and CBX2 ranked as the top three genes with genetic alterations, which altered in 25, 15, and 15% of the queried ESCA samples, respectively. The primary alteration type in these samples was the enhanced mRNA expression (Figure 5A). Besides, we calculated the correlations of CBXs with each other in GEPIA online tool for correlation analysis in ESCA patients and included Pearson’s correlation. The results showed that there were co-expression associations between the following CBX proteins: CBX1 positively with CBX2/3/4/5/6/8, CBX2 positively with CBX1/3/4/5/6/8, CBX3 positively with CBX1/2/5/8, CBX4 positively with CBX1/2/5/6/8, CBX5 positively with CBX1/2/3/4/6/7/8, CBX6 positively with CBX1/2/4/5/7/8, CBX7 positively with CBX5/6, and CBX8 positively with CBX1/2/3/4/5/6 (Figure 5B). Furthermore, CBX1/2/3/8 were positively co-expressed with TP53 (Figure 5C).

Using cBioPortal, we found 50 genes were most associated with each CBX family member. Some genes were positively associated with CBX family members, whereas others were negatively associated with the proteins. We used Metascape for Gene Ontology (GO) analysis of each CBX protein for biological processes, cellular components, and molecular functions. Moreover, functions of CBX family members and their 50 frequently altered neighbor genes were analyzed by Kyoto Encyclopedia of Genes and Genome (KEGG) and protein-protein interaction (PPI) enrichment analyses in Metascape. As shown in Figure 6A, biological processes such as GO: 0015931 (nucleobase-containing compound transport), GO:0006260 (DNA replication), GO:0006281 (DNA repair), GO: 0000226 (microtubule cytoskeleton organization), and GO: 0030029 (actin filament-based process) were prominently regulated by the CBX family members mutations in ESCA. Cellular components, including GO: 0098687 (chromosomal region), GO: 0000781 (chromosome, telomeric region), GO: 0090734 (site of DNA damage), GO: 0005635 (nuclear envelope) and GO: 0034399 (nuclear periphery) were remarkably associated with the CBX family members alterations. Besides, CBX family members mutations also significantly affected the molecular functions, such as GO: 0003779 (actin-binding), GO: 0003682 (chromatin binding), GO: 0032138 (single base insertion or deletion binding), GO: 0003697 (single-stranded DNA binding), and GO: 0015932 (nucleobase-containing compound transmembrane transporter activity).

In KEGG analysis, 11 pathways including ko:03430 (mismatch repair), ko:04512 (ECM-receptor interaction), ko:04910(Insulin signaling pathway), ko:04310 (Wnt signaling pathway), and ko:04514(Cell adhesion molecules (CAMs) were associated with the functions of CBX family members mutations in ESCA (Figure 6A). Figure 6B presented the PPI network associated with the genes; it was mainly related to cell cycle, gene mutation, DNA replication, and epigenetics.

Accumulating evidence indicated that immune cell infiltrations are closely related to tumor progression and clinical outcome. In this study, we also explored the correlation between CBX family members and immune cell infiltration in ESCA by utilizing the TIMER database (Supplementary Figure S4). There was a positive correlation between CBX1 expression and the infiltration of macrophages (Cor = 0.171, p = 2.20e-02), and a negative correlation between CBX1 expression and the infiltration of neutrophils (Cor = −0.15, p = 4.45e-02). CBX2 expression was negatively associated with the infiltration of CD8⁺ T cells (Cor = −0.184, p = 1.35e-02). CBX3 expression was negatively associated with the infiltration of dendritic cells (Cor = −0.249, p = 7.38e-04). Similarly, the expression of CBX5 was positively associated with the infiltration of macrophages (Cor = 0.199, p = 7.27e-03). There was a negative correlation between CBX6 expression and the infiltration of CD8⁺ T cells (Cor = -0.167, p = 2.51e-02), and a positive correlation between CBX6 expression and the infiltration of macrophages (Cor = 0.263, p = 3.62e-04). There was a positive correlation between CBX7 expression and the infiltration of B cells (Cor = 0.254, p = 5.87e-04), CD4⁺ T cells (Cor = 0.183, p = 1.43e-02) and macrophages (Cor = 0.283, p = 1.20e-04). The expression of CBX8 was positively associated with the infiltration of B cells (Cor = 0.222, p = 2.82e-03), and negatively associated with the infiltration of neutrophils (Cor = −0.164, p = 2.78e-02), and dendritic cells (Cor = −0.375, p = 2.09e-07).

TCGA ESCA project in level 3 HTSeq-FPKM format RNAseq data and clinical data showed that CBX1, CBX2, CBX3, CBX4, and CBX8 were negatively related to the infiltration of most immune cells, while CBX7 was positively associated with the infiltration of most immune cells (Figure 7). We also analyzed the relationship between CBX1/3/7/8 expression levels and immune cell infiltration in ESCA. Patients in the high CBX1 expression group (n = 173, data from TCGA database) presented a decrease in the numbers of infiltrating T cells, B cells, CD8 T cells, neutrophils, and Th17 cells (Supplementary Figure S5A). The group with higher expression of CBX3 showed a decrease in immune cell infiltration, including T cells, B cells, cytotoxic cells, CD8 T cells, neutrophils, and Th17 cells (Supplementary Figure S5B). The higher CBX7 expression group presented higher infiltration of immune which can kill tumor cells, including T cells, B cells, Cytotoxic cells, CD8 T cells, NK cells, iDC cells, Treg cells, and Th17 cells (Supplementary Figure S5C). Patients in the high CBX1 expression group indicated a decrease in the numbers of infiltrating T cells, cytotoxic cells, neutrophils, iDC cells, macrophages cells, and Treg cells (Supplementary Figure S5D).

These indicated that CBX family members may interact with immune cell infiltration, and then influence the outcomes of ESCA patients.

Aberrant epigenetic regulation has been reported to promote the pathogenesis and progression of ESCA. The methylation of lysine 27 on histone H3 (H3K27me3) is a chromatin marker associated with nucleosome condensation and silencing of gene expression (Adema and Colla, 2022). PCR1 and PCR2 played critical roles in establishing and maintaining the H3K27me3 mark. PRC2 comprised three core components (EZH2, SUZ12, and EED). EZH2 has hmtase activity, which is maintained by the presence of SUZ12 and EED (Erokhin et al., 2021). We used TCGA databases to analyze the expression levels of PRC2 components in ESCA. As shown in Supplementary Figure S6A, the transcriptional levels of EZH2, SUZ12, and EED were significantly elevated in ESCA tissues compared to normal tissues. We further used TCGA databases to analyze the association of CBX family members with PRC2 components in ESCA patients. The results showed that except CBX7, other CBX family members were positively correlated with PRC2 components, among which CBX1/3/8 had the highest correlation. (Supplementary Figure S6).

To validate the finding in the mentioned databases and further reveal which CBX members play a crucial role in the progression of ESCA, we used real-time PCR to detect all the mRNA expression of CBX family members using clinical samples. The analysis results showed that CBX1, CBX2, CBX3, CBX4, CBX5, CBX6, and CBX8 were significantly higher in ESCA tissues compared to normal esophageal tissues, which is consistent with the results of bioinformatics analysis mentioned previously (Figures 8A–F,H). However, CBX7 was downregulated in tumor tissues in half of collected ESCA patients, although there was no difference in all patients (Figure 8G).

We also investigated the protein expression of CBXs in ESCA tissues using western blotting and immunofluorescence staining. Western blotting results showed that the protein expression of CBX1, CBX3, and CBX8 were much higher in ESCA tissues than precancerous tissues (Figure 9). Immunofluorescence staining showed that the expression of CBX1, CBX3, and CBX8 were significantly higher in ESCA than adjacent normal tissues, which corroborated the results of western blot assay (Figure 10).