Mechanistic insights into super-enhancer-related genes as prognostic signatures in colon cancer

Data source

The GSE39582 dataset from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) database was defined as the training set for COAD, which contains 566 patients with COAD and 19 control population. The TCGA (Cancer Genome Atlas database, http://cancergenome.nih.gov) database was defined as the testing set for 467 COAD tissues and 41 normal tissues. The SERGs were derived from the SEdb database for HCT116 and HT29 cell lines and are shown in Supplementary Table 1.

Construction of WGCNA

Weighted Gene Co-expression Network Analysis (WGCNA) is performed as a systematic biology algorithm designed to construct gene co-expression networks and elucidate gene correlations across multiple dimensions. In this study, we utilized the R package “WGCNA” to estimate the COAD modules of correlated genes. Before the analysis, outliers were filtered out using the cutreeStatic function found within WGCNA. Subsequently, with a soft threshold power of 6 (b = 6), we generated the adjacency matrix to optimally fit the network structure. Pearson correlations were calculated between gene expression levels to construct the correlation matrix of genes, which established the connectivity between the nodes. Utilizing a hierarchical clustering dendrogram of this matrix, we built a topological overlap matrix to segregate different modules that follow similar gene expression patterns. The module eigengene (ME) expression profiles were then identified by amalgamating the expression profiles of each module, aiming to uncover a link between ME and clinical status. Candidate modules demonstrating significant correlation coefficients with clinical traits were subsequently shortlisted.

Enrichment analysis of SERGs

The genes from positive correlation coefficient modules of COAD in WGCNA and SERGs were intersected and visualized by using the Jvenn. Then, we analyzed functional annotations of these shared genes using GO and KEGG enrichment analysis.

Establish a prognostic risk model for SERGs

Univariate and multivariate Cox analysis, as the important means of determining prognostic-related SERGs in COAD, were performed for identification. Then, using the prognostic-related prediction formulas obtained by multivariate Cox regression analysis, the prognostic model was constructed in the R survival package. Afterwards, the high- and low-risk groups were differentiated by Kaplan Meier (K-M) survival analysis. Subsequently, the predictive value of the prognostic model was evaluated by the receiver operating characteristic area (ROC).

Immune profile analysis

We conducted a series of analyses to identify the differences between high- and low-SERGs risk groups in immune-cell infiltration. The abundance of 28 immune-cell types was determined using the ssGSEA. The Wilcoxon rank-sum test was applied to evaluate differences in immune cell proportions.

Single-cell analysis

We obtained single-cell RNA sequencing (scRNA-seq) data of COAD patients from the GEO database numbered GSE200997 and performed analysis using the ‘Seurat’ package. Subsequently, we identified immune cell clusters using the ‘FindNeighbors’ and ‘FindClusters’ functions. To evaluate the expression of key genes within the same type of immune cells between the tumor and normal groups, we utilized the ‘FindAllMarkers’ function.

Chemotherapeutic sensitivity

To predict the sensitivity of chemotherapeutic drugs between high- and low-SERGs risk groups, we utilized the pRRophetic R package to project drug sensitivity values (IC₅₀) by constructing a ridge regression model [13]. We selected several common anticancer drugs for this analysis, including camptothecin, docetaxel, gefitinib, gemcitabine, pazopanib, and sunitinib.

Hub genes validation in CC

To confirm the hub gene in CC, we performed the expression of ten genes on GSE39582 and TCGA-COAD via the Wilcoxon rank-sum test and ggplot2 with cutoff P-value < 0.05, respectively. The intersection of hub genes in the GSE39582 and TCGA-COAD was identified. The hub genes were identified using the GEPIA2 database based on “overall survival,” with the cutoff values determined using the “median” group. P < 0.05, P (HR) < 0.05 were significant. Subsequently, we acquired immunohistochemistry (IHC) images of CC and normal tissues via the Human Protein Atlas (HPA) portal [14].

Cell culture

CC cell lines (HCT116 and HT29 cells) and NCM460 cells were acquired from the Cell Bank of the Chinese Academy of Sciences. HCT116, HT29, and NCM460 cells were grown in the DMEM medium (10%-FBS) in a 37°C incubator with 5% CO₂.

qRT-PCR validation of the key gene

We isolated the total RNA from cells via TRIzol reagent (Invitrogen, USA). For cDNA synthesis, a cDNA Synthesis kit (Invitrogen, Thermo Fisher Scientific Inc., USA) was performed to reverse the transcription reaction into cDNA. The relative mRNA levels were assessed by the 2^−ΔΔCq calculation method and normalized by GAPDH mRNA expression. Primers were as follows: LZST2 forward primers: GGTGGCCCTATGACTTGG, reverse primers: AGCGGTGGGGAATGAAG. S100A11 forward primers: ATGGCAAAAATCTCCAGCCCT, reverse primers: TGTGAAGGCAGCTAGTTCTGTA. CYP2S1 forward primers: GCGCTGTATTCAGGGCTCAT, reverse primers: CTTCCAGCATCGCTACGGTT. PSMG1 forward primers: TCCTTTCCTGAGAGCCCTAAAA, reverse primers: TGTTCTAGCAATGGACAACACG. DCBLD2 forward primers: ATGTGGACACACTGTACTAGGC, reverse primers: CTGTTGGGATAGGTCTGTGGG. GAPDH forward primers: GGAAGCTTGTCATCAATGGAAATC, reverse primers: TGATGACCCTTTTGGCTCCC.

Western blot assay

CC cell lines (HCT116 and HT29 cells) and NCM460 cells were lysed in ice-cold RIPA lysis buffer (Beyotime Inc., China, P0013B). The quantification of protein content was achieved by utilizing the BCA Protein Assay Kit (Epizyme Biotech, China, ZJ101). The proteins were resolved via SDS-PAGE and subsequently transferred onto PVDF membranes. To block non-specific binding, the membranes were incubated with a solution containing 5% milk. Subsequently, the membranes were incubated with the corresponding primary antibodies overnight at 4°C. The horseradish Peroxidase secondary antibodies were applied and detected using an ECL solution (New Cell Molecular Biotech, China, P10200). Primary antibodies against LZTS2 (15677-1-AP) and GAPDH (60004-1-Ig) were sourced from Proteintech.

Statistical analysis

Statistical analyses were conducted using R (version 4.3.1) and GraphPad Prism. The T-test was performed to reveal the statistical differences between the two groups and One-way ANOVA analysis was for comparisons involving multiple groups. Statistical significance was determined at P < 0.05.

Data availability statement

The raw data are encompassed within the article and supplementary material.

Identification of SERGs

366 and 859 SERGs were identified in HCT116 and HT29 cell lines from the SEdb database, respectively. After removing the duplicates, a total of 974 SERGs were obtained.

Co-expression modules in LUAD and PAH

To identify the module genes associated with the disease, 23 modules were generated in GSE39582 by WGCNA, and different colors represented different modules. Then, we mapped the heat map, which could assess the association between modules and the disease based on the Spearman correlation coefficient (Figure 1). Considering the association of SEs with gene expression in tumors, most of the related genes may probably function as oncogenes [15]. Given that SEs enrich transcription factors to enhance gene expression, we specifically selected the modules “lightgreen”, “darkgrey”, “lightyellow”, “magenta”, “royablue”, “cyan”, “brown”, “darkturquoise”, “greenyellow”, and “blue” as COAD-related modules, which were positively correlated with COAD (lightgreen: r = 0.09, p = 0.03, genes = 209; darkgrey module: r = 0.25, p = 2e−09, genes = 147; lightyellow module: r = 0.25, p = 1e−09, genes = 203; magenta module: r = 0.28, p = 1e−11, genes = 654; royablue module: r = 0.25, p = 1e−09, genes = 191; cyan module: r = 0.25, p = 7e−10, genes = 398, brown module: r = 0.29, p = 7e−13, genes = 1612, darkturquoise module: r = 0.3, p = 2e−13, genes = 1770, greenyellow module: r = 0.15, p = 4e−04, genes = 996, blue module: r = 0.21, p = 4e−07, genes = 169).

Figure 1. Consensus module analysis of CC using WGCNA. (A) Scale-free topology model fit (R² > 0.98) and mean connectivity. (B) Association between the gene modules. (C) Cluster dendrogram for CC. (D) Heatmap analysis of modules and clinical features in CC.

Enrichment analyses of shared genes

A total of 406 genes at the intersection of WGCNA positive-correlated modules and SERGs were considered to be connected with the pathogenesis of COAD (Figure 2A). GO and KEGG analyses were utilized to identify 406 genes’ biological functions and essential pathways. The results indicated that the biological process (BP) was mainly enriched in tube morphogenesis and vasculature development, etc. The cellular component (CC) was primarily enriched in focal adhesion, etc. The molecular function (MF) was especially involved in cell adhesion molecule binding, etc., and KEGG analysis was primarily related to Pathways in cancer and PI3K-Akt signaling pathway (Figure 2B).

Figure 2. SERGs screening and function enrichment analysis. (A) Identification of pivotal SERGs in COAD. (B) GO and KEGG analysis of SERGs.

Establishment of the SERGs prognostic model

To further assess the SERGs associated with the prognosis of COAD and construct the model, we initially conducted univariate Cox analysis, considering genes with P < 0.01 (Supplementary Table 2). The P-value of the model (S100A11, LZTS2, RASA3, CYP2S1, ETFB, ZNF552, PSMG1, GJC1, NXN, and DCBLD2) was 2.3608e-09, with an AIC of 2225.98 by multivariate Cox analysis (Figure 3A). Even though the P-values for LZTS2, PSMG1, GJC1, NXN, and DCBLD2 exceed 0.05, it’s noteworthy that the model still exhibits the lowest AIC value among the tested models. We then proceeded to create a ten-gene SERGs model, which was formulated using the expression levels of each gene and their respective coefficients: risk score = (0.362746416 × S100A11) + (0.361912228 × LZTS2) + (0.495106924 × RASA3) + (−0.284872577 × CYP2S1) + (−0.274930212 × ETFB) + (−0.350131731 × ZNF552) + (−0.138081346 × PSMG1) + (0.395225421 × GJC1) + (−0.153361139 × NXN) + (0.203044412 × DCBLD2). Among the seven genes (S100A11, LZTS2, RASA3, GJC1 and DCBLD2) were classified as risk-related genes (HR > 1), while CYP2S1, ETFB, ZNF552, PSMG1 and NXN were protective genes (HR < 1). The risk score classified patients into high- and low-risk groups using the median as the threshold. To validate the signature of SERGs, we calculated the risk scores of patients in the TCGA-COAD dataset. The survival status in the two groups is shown in Figure 3B and gene expression in Figure 3C. Patients in the high-risk group suffered the worse OS. The validation results were largely consistent with those obtained from the GSE39582 dataset: the high-risk group endured worst OS (Figure 3D). Additionally, the AUCs at 1, 3, and 5 years for OS were 0.631, 0.684, and 0.681 in the GSE39582 dataset, respectively. The AUCs for OS at 1, 3, and 5 years were found to be 0.681, 0.702, and 0.633 in the TCGA-COAD dataset, respectively (Figure 3E).

Figure 3. Identification and validation of the ten-SERGs risk model. (A) The forest plot of ten-SERGs prognostic model. The risk score distribution, survival status, and heat map of ten SERGs in (B) training set and (C) testing set. (D) Patients in low-risk groups had longer OS in the training set and testing set. (E) The ROC analysis of the SERGs risk model in the training set and testing set.

Immune cell infiltration with SERGs risk group

We investigated the correlation between SERGs risk score and immune cell infiltration by ESTIMATE algorithms and ssGSEA. SERGs high-risk group had greater ESTIMATE score, Immune score, and Stromal score levels (P < 0.001) (Figure 4A). ssGSEA analysis revealed that the high-risk group had significantly lower levels of activated CD8 T cells, CD56 bright natural killer cells, CD56 dim natural killer cells, gamma delta T cells, memory B cells, monocyte, Type 17 T helper cells than those in the low group. However, the levels of effector memory CD8 T cells, Macrophages, Mast cells, natural killer T cells, natural killer cells, regulatory T cells, and Type 1 T helper cells were significantly higher than those in the low group (Figure 4B). The immune cells assessment between the two groups is shown in Figure 4C. Between the 28 immune cells, macrophage was positively correlated with myeloid-derived suppressor cells (r = 0.87) and regulatory T cells (r = 0.82) (Figure 4D).

Figure 4. Analysis of immune cell infiltration with SERGs risk group. (A) The ESTIMATE score, Immune score, and Stromal score between high- and low- SERGs risk groups were compared. (B) Violin diagram of 28 type immune cells in two groups. (C) Immune cells assessment between two groups. (D) Analysis of correlation in 28 type immune cells. (E) Correlation analysis of ten prognostic biomarkers (S100A11, LZTS2, CYP2S1, ZNF552, PSMG1, GJC1, NXN and DCBLD2) and immune cells. (F) Correlation analysis of risk score and immune cells.

Furthermore, the correlation analysis conducted on the ten biomarkers and the immune cells indicated a strong association between ZNF552 and S100A11 with immune cells (Figure 4E). The risk score was positively correlated with regulatory T cells, immature B cells and natural killer cells, and negatively correlated with memory B cells and Type 17 T helper cells (Figure 4F). These findings demonstrate the significant difference between the two groups in the immune microenvironment and are closely associated with ZNF552 and S100A11.

Single-cell analysis

The following major cell types were characterized: CD8 + T cells, CD4 + T cells, Plasma cells, epithelial cells, macrophages, B cells, Goblet cells, Natural killer cells, fibroblasts and endothelial cells (Figure 5A). S100A11 was mainly distributed in macrophages, endothelial cells, and fibroblasts. LZTS2 was mainly distributed in endothelial cells and fibroblasts (Figure 5B, 5C).

Figure 5. Single-cell analysis. (A) 10 types of cells were clustered. (B) The ten identified SERGs markers expressions were identified in CC single-cell clusters. (C) A bubble plot was employed to visually represent the gene expression characteristics.

SERGs were predictive to chemotherapy

We employed the pRRophetic algorithm to estimate IC₅₀ values, enabling the prediction of distinct chemotherapy responses between high- and low-SERGs risk groups. Based on the GSE39582 database, we found that the high-SERGs risk group had a low IC₅₀ of the anticancer drugs compared to high-SERGs risk group such as docetaxel (P < 0.02), gefitinib (P < 0.001), pazopanib (P < 0.001), sunitinib (P < 0.001), which means that the above drugs were more sensitive to the high-SERGs risk group (Figure 6).

Figure 6. The pRRophetic algorithm predicted the IC₅₀ values for six anti-cancer drugs.

Identification of key gene expression in COAD

To further identify the key gene, we analyzed the 10 genes between cases and controls on GSE39582 and validation on the TCGA-COAD dataset. Eight genes were selected, including S100A11, LZTS2, CYP2S1, ZNF552, PSMG1, GJC1, NXN and DCBLD2 (Figure 7A). Meantime, S100A11, LZTS2, CYP2S1, PSMG1 and DCBLD2 were found on validation cohorts TCGA-COAD dataset (Figure 7B). Finally, we found that the five genes (S100A11, LZTS2, CYP2S1, PSMG1 and DCBLD2) were up-regulated in both GSE39582 and TCGA-COAD datasets with P < 0.05. Then, we performed qRT-PCR and found that LZTS2 was significantly up-regulated in HCT116 and HT29 cell lines compared to normal NCM460 cells (Figure 7C). The expression of S100A11, CYP2S1, PSMG1, and DCBLD2 were inconsistent with predicted results. Specifically, the five genes were queried in the GEPIA2 database, and LZST2 was selected based on “overall survival.” (Supplementary Figure 1). Notably, the WB experiment showed that the levels of LZTS2 decreased in HCT116 and HT29 cell lines than in NCM460 cells (Supplementary Figure 2). The high LZTS2 expression was related to an unfavorable OS by prognostic analysis (Figure 7D). Finally, the immunohistochemical findings retrieved from the HPA database demonstrated elevated LZTS2 protein expression within the CC tissue (Figure 7E).

Figure 7. Identification of LSZT2 mRNA expression. (A) The expressions of ten genes in discovery cohorts. (B) The expression of ten genes in validation cohorts. (C) The LSZT2 expression level in NCM460 cells and CC cell lines. (D) The prognostic analysis of LSZT2 in CC patients. (E) IHC staining of LSZT2 in normal and CC tissues from the HPA database.