Unveiling the prognostic significance of SOX5 in esophageal squamous cell carcinoma: a comprehensive bioinformatic and experimental analysis

GEO datasets

GEO (http://www.ncbi.nlm.nih.gov/geo) contains microarrays, chips and gene expression data [19]. GSE23400 (53 ESCC samples and 53 matched normal samples) was used in this study. All downloaded data were matched with clinical information, standardized and log2 translated.

TCGA datasets

The RNA-seq data of SOX5 and the clinical information of TCGA-ESCC were from the UCSC Xena browser (https://xenabrowser.net). The samples data including the information of age, gender, clinical stage, lymph node metastasis, neoplasm histologic grade and overall survival (OS) were analyzed. Furthermore, the samples were divided into two groups in line with median SOX5 expression, Gene Set Enrichment Analysis (GSEA) was conducted. At the same time, the sample contains RNA-seq data of SOX5 and methylation data were filtered. The Champ package was used to pre-process the methylation data, and the CpG sites in the SOX5 promoter region (TSS200, TSS1500) were selected. The correlation between the beta value of these CpG sites and SOX5 expression was calculated using Pearson analysis. CpG sites with SNP were ignored.

Identification of DEGs

DEGs between ESCC and normal samples were identified by Limma package (|logFC|>1 and p-value<0.05). First, we preprocessed the data to remove biases and normalize the expression values. Next, we constructed a design matrix to model the experimental design and define the comparisons of interest. Using this design matrix, a linear model was fitted to the data, considering the experimental factors and estimating treatment effects. We then defined contrasts to specify the desired comparisons between conditions. Statistical testing was performed using moderated t-statistics to evaluate the significance of differential expression. Multiple testing correction was applied, and DEGs were identified based on adjusted p-values and fold change thresholds. The resulting DEG list was annotated with relevant information, such as gene symbols and pathways. Finally, visualizations and analyses were generated to aid in the interpretation of the DEGs. Based on Pearson correlation analysis, the correlation of DEGs and SOX5 expression was analyzed.

ChIP-seq data analysis

The results of bed peaks of ChIP-seq data for SOX5 were provided by Cistrome DB (http://cistrome.org/db/#/). The downloaded narrow peak data were annotated and filtered with a threshold value of Q-value and IDR < 0.01. The peaks within the range of 1 kb upstream or downstream of transcription start site (TSS) were retained, and the genes on these peaks were defined as target genes. Dreme was used to identify significant motif sequences in peak sequences.

Enrichment analysis

Biological functions of target genes were conducted by clusterProfiler package in Bioconductor. Gene Ontology (GO), Kyoto Encyclopedia of Genes and Genomes (KEGG), and Gene Set Enrichment Analysis (GSEA) were performed to identify enriched biological processes, molecular functions, cellular components, and pathways associated with a set of genes of interest. GO enrichment analysis compared the genes against a background set using statistical tests and adjusted p-values or false discovery rate (FDR) correction. Similarly, KEGG enrichment analysis compared genes to KEGG pathways, while GSEA ranked genes and evaluated enrichment scores through permutation testing. Enriched terms and pathways were identified based on predefined thresholds, providing insights into gene function and pathway associations.

Protein-protein interaction (PPI) network construction

STRING (http://string-db.org, version 11.0) database was used to construct the PPI network of target genes. Next, we use Cytoscape software to visualize the protein-protein interaction network and describe the relationship between different proteins by line segments, in which the more line segments represented the more significant the interaction between proteins.

Hub genes selection and analysis

Nodes were ranked in the PPI network with respect to their network features through Cytoscape, which is the cytoHubba plug-in [20]. Here, the maximal clique centrality (MCC), and edge percolated component (EPC) were applied as the criteria for ranking the top 5 hub genes. The overlapping top 5 genes in these two topological methods were regarded as hub genes for further analysis [21].

Immune infiltration analysis

The 22 immune cell types in samples were assessed using the CIBERSORT algorithm with the criteria of p-value <0.05. The proportion of each kind of immune cell in the samples was calculated. The Wilcoxon rank-sum test was utilized to determine the percentage of immune-infiltrating cells between ESCC and control samples. After inclusion of all significant cells into the Cox model, co-expressed heatmap was drawn by combining with Pearson correlation analysis, and the correlation of various immune cells and hub genes was displayed.

Cell lines and cell culture

Human ESCC cell lines KYSE150, KYSE30, and normal esophageal epithelial cell lines (Het1A) were cultured in RPMI-1640 with 10% FBS (GIBCO, Waltham, MA, USA) at 37° C under 5% CO₂.

Plasmid construction and cell transfection

The full-length of SOX5 cDNA was generated and inserted into the pAd-CMV-MC5 to construct SOX5 overexpression vector (Ad-SOX5) [22]. Lipofectamine 2000 (Thermo Fisher Scientific, Cleveland, OH, USA) was applied for cell transfection.

Western blot

After lysing, protein of cells was separated by SDS-PAGE and transferred onto PVDF membrane (Millipore, USA). Primary antibodies were incubated on the membrane at 4° C overnight: SOX5 (1:1000, Abcam), PCNA (1:1000, Abcam), RRM2 (1:1000, Abcam), AURKB (1:1000, Abcam), MCM4 (1:2000, Abcam), MCM7 (1:1000, Abcam), GAPDH (1:1000, Abcam). The results were determined with a chemiluminescence kit (Thermo Fisher Scientific, Cleveland, OH, USA).

Cell proliferation assay

For CCK-8 assay, treated cells in 96-well culture plates (1.5 × 10³ cells/well) were treated with 10 μL CCK-8 solution (Dojindo, Kumamoto, Japan) at different time points. After 2 h, the optical density (OD) at 450 nm was detected.

For EdU assay, treated cells in 96-well plates (1.5 × 10⁴ cells/well) were incubated with EdU working solution (Cat.C10310-1, Ruibo Biotech., Guangzhou, China) for 2 h. After fixation, permeabilization, cells were stained with 1x Apollo solution. Results were obtained under fluorescence microscope.

Flow cytometry

Annexin V-FITC apoptosis detection kit (Oncogene Research Products, Boston, MA, USA) was used to determine the apoptosis of cells.

Statistical analysis

R V4.0.2 was used to perform statistical analysis. Survival rates were determined by Kaplan-Meier survival curves. Pearson’s correlation test was utilized to detect the extent of correlation between the expression of SOX5 and the methylation/DEGs, and hub genes. The differences between two groups were analyzed by GraphPad Prism 8. p< 0.05 indicates a statistically significant difference.

Low expression of SOX5 in ESCC is associated to survival

To explore the SOX5 expression in human ESCC tissues, we downloaded microarray gene profiling datasets GSE23400 (53 ESCC samples and 53 matched normal samples) from the GEO. We analyzed SOX5 expression in datasets, and the results are shown in Figure 1A, indicating that SOX5 was downregulated in ESCC tumor tissue compared with normal tissue (p=0.028). The results indicated that SOX5 might play an important role in ESCC.

Figure 1. Low expression of SOX5 in ESCC is associated to survival. (A) Expression boxplots of SOX5 in ESCC tumor tissue and adjacent normal tissue from GSE23400. (B) Survival analysis of SOX5 in ESCC samples from TCGA-ESCC. (C) Univariate and multivariate cox regression analyses of clinicopathological characteristics for OS in ESCC samples from TCGA-ESCC. HR, hazard ratio; CI, confidence interval.

Therefore, we screened 80 ESCC samples in TCGA-ESCC from the XENA. First, the relationship between SOX5 expression and survival outcomes in the ESCC cohort was evaluated (Figure 1B, p=0.0076). Furthermore, SOX5 low expression group had a shorter survival time than the high expression group. Additionally, univariate analysis and multivariate analyses of Cox regression model were performed. In univariate analysis, low SOX5 expression and gender were significantly correlated to OS of ESCC patients (Figure 1B). In multivariate analysis, low SOX5 expression was significantly correlated to OS in patients with ESCC (Figure 1C). The results of COX regression are shown in Supplementary Figure 2.

Results of GSEA in ESCC

To gain a deeper understanding of the biological functions of SOX5 in ESCC, we conducted GSEA. Similarly, the samples from TCGA-ESCC cohort were divided into two groups in line with the median SOX5 expressions in the cohort as a threshold. As the results showed, low expression of SOX5 group was mainly enriched in the pathway of regulated necrosis, the NLRP3 inflammasome, formation of the cornified envelope and PD 1 signaling (Figure 2A–2D).

Figure 2. Results of GSEA in ESCC. (A–D) GSEA enrichment analysis of SOX5 in TCGA-ESCC.

SOX5 expression is negatively correlated with its methylation in ESCC

To understand the potential regulatory effect of methylation on SOX5 expression, we downloaded the TCGA-ESCC methylation chip data from XENA. After filtering, 74 ESCC samples with both RNA-Seq data of SOX5 and methylation data were obtained. The CpG sites on the SOX5 gene were filtered, the CpG sites located in the promoter region were retained, and 57 CpG sites were obtained. The association between SOX5 expression and the β-value of CpG sites was identified using the Pearson’s correlation analysis. Results (Figure 3A, 3B) indicated a negative correlation between SOX5 expression and 2 sites: cg01198491 (Pearson R=-0.37; p<0.01), and cg06141624 (Pearson R= -0.3; p<0.05).

Figure 3. SOX5 expression was negatively correlated with its methylation in ESCC. The correlation of SOX5 expression and beta-value of cg01198491 (A) and cg06141624 (B).

Identification of DEGs in ESCC

Via Limma, DEGs in the GSE23400 dataset were identified (Figure 4A, 4B). By screening, 187 upregulated and 207 downregulated DEGs were obtained. The results were presented with the volcano map and heatmap. The correlation between the above DEGs and SOX5 was calculated by Pearson correlation coefficient with a PCC > 0.2 and a p-value < 0.05. Figure 4C shows that 122 genes were significantly correlated with SOX5.

Figure 4. Identification of DEGs in ESCC. (A) Volcano map and (B) heatmap of DEGs between ESCC samples and control samples from GSE23400. (C) Correlations between SOX5 and DEGs in ESCC were analyzed using Pearson-correlation analysis.

ChIP-seq data analysis of target genes of SOX5

To further explore the correlated genes of SOX5, we predicted target genes of SOX5 using ChIP-seq data assay and predicted the SOX5 binding sites. Distribution of binding peaks in the SOX5 ChIP-seq experiment from Cistrome DB were analyzed. The peaks in the range of I kb of upstream and downstream of TSS were reserved. Approximately 6324 peaks were located on 4951 genes (Figure 5A–5C). We took the intersection of these genes and 122 significantly correlated DEGs, and finally 28 target genes were obtained (Figure 5D). And DREME (https://meme-suite.org/meme/tools/dreme) online software was used to predicted SOX5 binding motifs (Figure 5E).

Figure 5. ChIP-seq analysis of target genes of SOX5. (A) Heatmap of ChIP binding to TSS regions. (B) Average tag density profiles. (C) Pie plot of genomic regions of peaks. (D) The intersection of the results of ChIP-seq data analysis and significantly correlated DEGs. (E) Top5 predicted motifs by DREME.

Gene ontology and pathway functional enrichment analysis

A total of 28 target genes were conducted with GO and KEGG enrichment analysis. Genes were enriched in DNA replication, G1/S transition of mitotic cell cycle, and cell cycle G1/S phase transition in the biological process (BP) analysis (Figure 6A). For molecular function (MF) analysis, genes were primarily linked to the ATPase activity, helicase activity, and catalytic activity, acting on DNA (Figure 6B). The enrichment analysis of cell component (CC) indicated that genes were related to spindle, chromosomal region, and mitotic spindle (Figure 6C). For KEGG enrichment analysis, genes were enriched mainly in DNA replication, cell cycle and ferroptosis (Figure 6D). These analyses suggested that target genes had a crucial effect on ESCC.

Figure 6. Gene ontology and pathway functional enrichment analysis. (A) BP analysis. (B) MF analysis. (C) CC analysis. (D) KEGG pathway. Functional and pathway enrichment were presented in bubble charts, ranking the top 10 according to adjusted P value.

Hub genes identification

The potential relationships among 28 target genes were detected by STRING database, and PPI network was established (minimum required interaction score > 0.4) (Figure 7A). Next, two topological analysis methods, MCC and EPC, were used to rank the top 5 genes in the PPI network (Figure 7B). The results showed that 5 overlapping hub genes (PCNA, RRM2, AURKB, MCM4 and MCM7) were identified. And ChIP-seq profiles of SOX5 at promoter region of hub genes were presented in Figure 7C.

Figure 7. Hub genes identification. (A) The PPI network of DEGs was established by STRING. (B) The top 5 nodes of PPI network evaluated in cytoHubba MCC and EPC arithmetic. (C) ChIP-seq profiles of SOX5 at promoter region of PCNA, RRM2, AURKB, MCM4, MCM7.

Pathway enrichment analysis of hub genes

The expressions of 5 hub genes were significantly higher in ESCC samples compared to normal samples (Supplementary Figure 1). The results of pathway analysis (Figure 8A) showed that 3 hub genes (PCNA, MCM4 and MCM7) were related to DNA replication and cell cycle pathway (Figure 8B, 8C). These findings are consistent with the molecular pathways implicated in ESCCC carcinogenesis. Further, the previous study indicated that the elevated PCNA+ tumor-associated macrophages in early recurrence of BC [23]. MCM7 has been found to be involved in the immune response of the spleen, thereby influencing tumor development [24]. These results implied that hub genes involve in the tumor progression and immune response.

Figure 8. Pathway enrichment analysis of hub genes. (A) KEGG pathway analysis. KEGG pathway annotations of DNA replication signaling (B) and cell cycle signaling pathways (C).

Immune infiltration analysis

Therefore, we investigated whether hub genes were correlated with immune infiltration levels in ESCC. The composition of 22 immune cells (Figure 9A) and the distribution of immune cells (Figure 9B) in samples uncovered an important relationship between ESCC and immune cell infiltration. In addition, Figure 9C shows that there is a noteworthy positive correlation between 5 hub genes and the infiltration of Macrophages.M0 and Macrophages.M1.

Figure 9. Immune infiltration analysis. (A) The composition of immune cells estimated in ESCC. (B) Differential expression of different types of immune cells between ESCC and normal samples. (C) The relationships among immune cells and between hub genes and immune cells.

Upregulation of SOX5 inhibited proliferation and promoted apoptosis of ESCC cells in vitro

We sought to elucidate the expression of SOX5 and 5 hub genes through studies on ESCC cells in vitro. Using Western blot, SOX5 was downregulated and 5 hub genes were upregulated in different ESCC cell lines, compared to Het1A cells (Figure 10A). Furthermore, SOX5 in ESCC cell lines was overexpressed (Figure 10B). Subsequently, upregulation of SOX5 showed inhibitory effects on ESCC cell proliferation (Figure 10C, 10D). Flow cytometric analysis and Western blot showed Upregulation of SOX5 elicited apoptosis of KYSE30 and KYSE150 cells (Figure 10E, 10F) (*p<0.05, **p<0.01, ***p<0.001).

Figure 10. Upregulation of SOX5 inhibited proliferation and increased apoptosis of ESCC cells in vitro. (A) Western blot assay of SOX5 and hub genes expression in ESCC and Het1A cells. (B) Western blot assay of SOX5 protein in control and overexpression ESCC cells. (C, D) CCK-8 assay (C) and EdU assay (D) results of different treated cells. (E) Western blot assay of Bcl2, Bax, cleaved-caspase 3 and caspase 3 protein in different treated cells. (F) Flow cytometry analysis of apoptosis of different treated cells. *P < 0.05, **P < 0.01.