A cohesin-associated gene score may predict immune checkpoint blockade in hepatocellular carcinoma.

Stromal antigen 1 (STAG1), a component of cohesion, is overexpressed in various cancers, but it is unclear whether it has a role in the transcriptional regulation of hepatocellular carcinoma (HCC). To test this hypothesis, here, we screened global HCC datasets and performed multiscale embedded gene co-expression network analysis to identify the potential functional modules of differentially expressed STAG1 co-expressed genes. The putative transcriptional targets of STAG1 were identified using chromatin immunoprecipitation followed by high-throughput DNA sequencing. The cohesin-associated gene score (CAGS) was quantified using the The Cancer Genome Atlas HCC cohort and single-sample gene set enrichment analysis. Distinct cohesin-associated gene patterns were identified by calculating the euclidean distance of each patient. We assessed the potential ability of the CAGS in predicting immune checkpoint blockade (ICB) treatment response using IMvigor210 and GSE78220 cohorts. STAG1 was upregulated in 3313 HCC tissue samples compared with 2692 normal liver tissue samples (standard mean difference = 0.54). A total of three cohesin-associated gene patterns were identified, where cluster 2 had a high TP53 mutated rate and a poor survival outcome. Low CAGS predicted a significant survival advantage but presaged poor immunotherapy response. Differentially expressed STAG1 co-expression genes were enriched in the mitotic cell cycle, lymphocyte activation, and blood vessel development. PDS5A and PDGFRA were predicted as the downstream transcriptional targets of STAG1. In summary, STAG1 is significantly upregulated in global HCC tissue samples and may participate in blood vessel development and the mitotic cell cycle. A cohesin-associated gene scoring system may have potential to predict the ICB response.

cohesin‐associated gene

cohesin‐associated genes scores

chromatin immunoprecipitation sequencing

cell type identification by estimating relative subsets of RNA transcripts

complete remission

estimation of STromal and immune cells in MAlignant tumor tissue using expression data

hepatocellular carcinoma

immune checkpoint blockade

multiscale embedded gene co‐expression network analysis

progressive disease

partial remission

stable disease

standard mean difference

single‐sample gene set enrichment analysis

stromal antigen 1

weighted gene co‐expression network analysis

Hepatocellular carcinoma (HCC), one of the most pernicious malignancies in the digestive system, comprises approximately 80% of liver cancer cases and poses a significant threat to the global population [1]. Risk factors for HCC vary from region to region, and nonalcoholic fatty liver disease has become the fastest‐growing cause of HCC in the United States (US) and the United Kingdom [2]. Although surgical resection, radiofrequency ablation, and transcatheter arterial chemoembolization have been used to treat HCC patients, there is a dearth of effective therapeutic strategies for progressed HCC [3, 4, 5]. Immunotherapy exhibits impressive prospects in treating HCC patients [6]; however, significant heterogeneity exists when responding to immune checkpoint blockade (ICB) treatment [7, 8]. Therefore, there is an urgent need to identify more effective indicators for predicting the response to ICB treatment [9]. Moreover, the ultimate pathogenesis of HCC must be identified, and more effort must be made to treat recurrent and metastatic HCC patients [10, 11].

The cohesin core subunit, stromal antigen 1 (STAG1), is encoded by the STAG1 gene, which belongs to the sister chromatid cohesion protein 3 family, and is ubiquitously expressed in the nucleus [12]. It is known that the cohesin protein complex is indispensable for sister chromatid cohesion; it also plays important roles in transcriptional regulation in addition to chromosome maintenance [13]. Interestingly, STAG1 exerts its function by binding to cohesin with the aid of the RAD21 cohesin complex component, thus generating a platform for the binding of other regulatory cohesin subunits [14]. It has been reported that cohesin can promote DNA damage repair [15], and cohesin mutations are commonly detected in cancers [16, 17, 18]. Moreover, inactivated STAG1 was found to attenuate the proliferation of bladder cancer and sarcoma cells [19]. Nonetheless, the role of the STAG1 transcriptional factor (TF) in HCC remains unknown.

Given the lack of research on STAG1, this study focused on identifying novel cohesin‐associated HCC phenotypes and exploring the underlying transcriptional regulatory mechanism, thus providing avenues for treating HCC patients.

Materials and methods

Cohesin‐associated gene signatures in the cancer genome atlas liver HCC cohort

The design route of this research study is displayed in Fig. S1. Level three fragments per kilobase of transcript per million fragments mapped (FPKM) data were downloaded from The Cancer Genome Atlas (TCGA) and were transformed into transcripts per million (TPM) values. As a pivotal component of cohesin, STAG1 exerts its function with other cohesin‐associated genes (CAGs). Therefore, the expression data of 17 CAGs were abstracted, including five subunits (RAD21, SMC1A, SMC3, STAG1, and STAG2), six regulatory factors (CDCA5, MAU2, NIPBL, PDS5A, PDS5B, and WAPL), and six regulation controlling factors (AURKB, CDK1, ESPL1, PLK1, PTPA, and SGO1) [20]. Single‐sample gene set enrichment analysis (ssGSEA) was used to quantify the cohesin‐associated signatures in each of the HCC patients, which were defined as the CAG scores (CAGS).

Identification of the CAG patterns in HCC patients

Principal component analysis (PCA) was conducted to determine the ability of the CAGs to discriminate between HCC and the adjacent normal liver tissue samples. HCC tissue samples were assigned to low, moderate, and high CAGS clusters by calculating their Euclidean distance based on CAGS [21].

Clinical prognosis and molecular characterization of different CAG patterns in HCC patients

HCC patients were assigned to either a high CAGS group or a low CAGS group, based on the median CAGS value. Kaplan–Meier survival analysis was used to compare the overall survival (OS) period of the high CAGS and low CAGS HCC patients and that of three CAG patterns. The mutation annotation format (MAF) files of the HCC patients were acquired from TCGA, and the somatic mutation landscapes of the three CAG patterns were compared. The immune scores, stromal scores, and tumor purity of the HCC patients were calculated using the Estimation of STromal and Immune cells in MAlignant Tumor tissue using the Expression data (ESTIMATE) algorithm [22, 23]. The infiltration levels of 22 immune cells were quantified in each HCC patient using the Cell type Identification by Estimating Relative Subsets of RNA Transcripts (CIBERSORT) algorithm [24, 25, 26, 27, 28]. The tumor microenvironment (TME) of different CAG patterns was compared using either the Wilcoxon test or the ANOVA. The metabolism‐related signatures in each HCC patient were quantified by PCA, and the metabolic discrepancies of the three CAG patterns were compared in a heatmap.

Prediction of ICB treatment response

To detect the potential ability of the CAGS in predicting ICB treatment response, we downloaded the expression data and clinical features of two ICB cohorts (i.e., the IMvigor210 urothelial cancer [29] cohort and the GSE78220 melanoma [30] cohort). The IMvigor210 and GSE78220 expression matrices were transformed into TPM data. Two patients in the GSE78220 cohort were excluded because they lacked survival data or they had previously received any other on‐therapy. For the duplicated probes in each matrix, only the mean expression values were preserved. The CAGS of each patient was quantified and was used to predict the OS condition and the ICB therapeutic response of the patients, which included progress disease (PD), stable disease (SD), partial remission (PR), and complete remission (CR).

In‐house immunohistochemistry (IHC)

HCC tissue microarrays (TMA) LVC1505 and LVC1531 were purchased from Fanpu Biotechnology Co., Ltd, Guilin, China (http://www.fanpu.com/). All the tumor specimens were pathologically diagnosed with primary HCC, and the patients did not receive drug interventions before surgery. IHC staining was performed using the anti‐STAG1 antibody (Abcam Co., Ltd, ab246988) according to the procedure previously reported [31, 32, 33, 34]. The dilution of the anti‐STAG1 antibody was 1 : 100. The patients provided written informed consent for the use of their tissue samples. This study has been approved by the Ethics Committee of The First Affiliated Hospital of Guangxi Medical University [No.2022‐KY‐E‐017]. The protocols conformed to the guidelines set by the Declaration of Helsinki.

Global HCC microarrays and RNA sequencing (RNA‐seq) datasets

Public HCC datasets, including microarrays and RNA‐seq matrices, were searched in the Gene Expression Omnibus (GEO) and ArrayExpress databases. The query keywords were as follows: HCC OR hepatocellular carcinoma OR liver cancer OR liver tumor. The inclusion criteria were as follows: (a) The specimens should be primary HCC tissue samples from humans, (b) the patients did not receive preoperative treatment, and (c) the sample size of each platform should be ≥ 6. TCGA was combined with the Genotype‐Tissue Expression (GTEx) project. We used the combat function to remove the batch effect of the same GEO platform and that of the integrated TCGA‐GTEx dataset. Data normalization was performed with the aid of the limma‐voom (The Walter and Eliza Hall Institute of Medical Research, Melbourne, Australia) package in r v3.6.1.

mRNA expression statuses of CAG and HCC differentially expressed genes (DEGs) based on global HCC data

The mRNA expression level of STAG1 and 16 other CAGs was calculated using the enrolled global HCC datasets. Standard mean difference (SMD) was chosen as a comprehensive index for appraising the relative expression status of CAG, where an SMD score > 0 or < 0 (P value < 0.05) represented an upregulated or downregulated trend, respectively. Kaplan–Meier OS analysis was performed to evaluate the prognostic value of STAG1 mRNA in the TCGA‐HCC cohort.

STAG1 co‐expression genes (CEGs)

STAG1 CEGs were identified from the processed global HCC datasets by calculating Pearson correlation coefficients, where coefficients ≤ −0.3 or ≥ 0.3 and P values < 0.05 represented negatively correlated STAG1 CEGs or positively correlated STAG1 CEGs, respectively. Weighted gene co‐expression network analysis (WGCNA) was conducted to determine the high‐CAGS HCC phenotype‐related module genes.

Multiscale embedded gene co‐expression network analysis (MEGENA)

We acquired the HCC differentially expressed STAG1 CEGs by intersecting the upregulated DEGs with the positively correlated STAG1 CEGs, and by overlapping the downregulated DEGs with the negatively correlated STAG1 CEGs. The functional networks of the HCC differentially expressed STAG1 CEGs were identified by analyzing their topological structure. Metascape was utilized to annotate the biological functions of the HCC differentially expressed STAG1 CEGs.

Putative transcriptional targets of STAG1 based on global chromatin immunoprecipitation sequencing (ChIP‐seq) data

A total of 58 Cistrome ChIP‐seq datasets were utilized to identify the candidate transcriptional target genes of the STAG1 factor. The cutoff value for filtering the putative targets of STAG1 was set at a score ≥ 3.0. The transcriptional regulatory relationships were further validated by peak annotations and Pearson co‐expression analysis.

Potential therapeutic agents for HCC

The CellMiner database is a repository of molecular and pharmacological data based on 60 cell lines from the National Cancer Institute (NCI‐60), which contains the transcriptome data of 60 cancer cell lines and over 100,000 natural products and chemical compounds [35]. Herein, both US Food and Drug Administration (FDA) approved drugs, and clinical trial drugs were selected for in‐depth analysis. To screen for putative anticancer therapeutic agents, we first analyzed the discrepancies of half maximal inhibitory concentration (IC50) in the STAG1 CEG high and low expression groups by targeting the STAG1‐enriched molecular pathways. We computed Pearson correlation coefficients between the IC50 values of the screened agents and the STAG1 CEG expressions. Computer‐aided drug‐protein molecular docking was conducted to fundamentally validate the potential therapeutic agents for treating HCC.

Results

Distinct CAG patterns in HCC patients with different clinical phenotypes

In this study, the CAGs showed a highly effective performance in differentiating HCC tissue samples from adjacent normal liver tissue samples (Fig. 1A). A total of three CAG patterns were identified from the HCC patients (Fig. 1B), where cluster 2 showed a worse OS outcome than any other clusters (Fig. 1C). Figure 1D summarizes the clinical characteristics of the different HCC clusters. We also quantified the CAGS of each HCC patient, and noticed that the low CAGS HCC patients displayed a significant survival advantage (Fig. 1E). Among the 17 CAGs, a total of 13 signatures were significantly overexpressed in global HCC tissue samples (Fig. 1F).

Fig. 1

STAG1‐mediated distinct prognostic and mutation phenotypes in HCC. Given that STAG1 is a pivotal component of cohesin, it and 16 other CAGs were collected to explore their biological functions. (A) Using STAG1 and 16 other CAGs, HCC tissue samples could be preferably differentiated from adjacent normal liver tissue samples. (B) HCC patients were assigned to three clusters based on their cohesin‐associated signatures. (C) HCC patients classified into cluster 2 were found to have a worse survival outcome than the patients in the other clusters. A log‐rank test was used to identify the survival discrepancy. (D) The clinical characteristics of different HCC clusters. (E) The CAG signatures in each patient were quantified using single‐sample gene set enrichment analysis. The calculated values, termed CAGS, were used to classify the HCC patients. The low CAGS HCC patients displayed a significant survival advantage. A log‐rank test was used to identify the survival discrepancy. (F) Dysregulated expression levels of STAG1 and 16 other CAGs in the global HCC data (sample size: T  = 3831, N  = 3102; T  = 3412, N  = 2653; T  = 2347, N  = 1832; T  = 3926, N  = 3428; T  = 2141, N  = 1706; T  = 3334, N  = 2729; T  = 3081, N  = 2481; T  = 3055, N  = 2451; T  = 3414, N  = 3036; T  = 384, N  = 305; T  = 3395, N  = 3023; T  = 397, N  = 296; T  = 3386, N  = 3008; T  = 3306, N  = 2701; T  = 3313, N  = 2692; T  = 3414, N  = 3036; T  = 425, N  = 340, where T and N stand for HCC and non‐HCC group, respectively). Distinct mutation phenotypes of the three HCC clusters are shown in G–I, where TP53 is identified as the most frequently mutated gene in the HCC patients of cluster 2. CAG, cohesin‐associated gene; CAGS, cohesin‐associated gene scores; HCC, hepatocellular carcinoma; STAG1, stromal antigen 1. ****P value < 0.0001; *P value < 0.05; ns, not significant.

Metabolic characterization and genetic alterations and TME landscapes of three CAG patterns

Given the essence of tumor metabolic reprograming in cancer progression [36], we first compared the metabolic characterization based on the identified CAG clusters. Notably, the HCC patients in CAG cluster 2, with a poorer prognosis than any of the other clusters, displayed prominently lower metabolism levels of drug, retinol, and xenobiotics (Fig. S2A–C). Additionally, distinct gene mutation phenotypes were identified among the three HCC CAG clusters (Fig. 1G–I), where TP53 was the most frequently mutated gene in cluster 2. Moreover, the mutation co‐occurrence between TP53 and other genes, such as TTN, MUC16, and CTNNB1, was more frequent in HCC CAG cluster 2 (Fig. S2D–F) than in the other clusters. As is shown in Fig. 2A, the STAG1 mRNA expression level was compared between different CAG clusters. Pertaining to TME, HCC CAG cluster 2 exhibited significantly lower immune and stromal scores (Fig. 2B,C) but possessed higher tumor purity (Fig. 2D). Additionally, HCC CAG cluster 2 showed higher infiltration levels of resting dendritic and T follicular helper cells (Fig. 2E,F) than the other clusters.

Fig. 2

Distinct TME characterization of HCC patients. According to the identified CAG clusters, the TME characteristics of different HCC patients were compared. A Wilcoxon test was used to compare the intercluster immune and stromal scores (sample size: Ncluster1 = 251; Ncluster2 = 85; Ncluster3 = 35). The STAG1 mRNA expression level was significantly higher in CAG cluster 2 (A). HCC patients in CAG cluster 2 with poor prognosis exhibited prominently lower immune scores and stromal scores (B–C), and possessed higher tumor purity (D) than the other HCC patients. An ANOVA method was used to compare the immune infiltration levels among three clusters (sample size: Ncluster1 = 251; Ncluster2 = 85; Ncluster3 = 35). HCC patients in CAG cluster 2 also showed higher infiltration levels of resting dendritic and T follicular helper cells (E–F). HCC, hepatocellular carcinoma; CAG, cohesin‐associated gene; TME, tumor microenvironment. ****P value < 0.0001; ***P value < 0.001; **P value < 0.01; *P value < 0.05; ns, not significant.

Potential ICB immunotherapy response prediction ability of CAGS

As is shown in Fig. S3, higher STAG1 expression level predicted better OS outcomes in patients who received immunotherapy. In light of the role of CAGS in predicting OS for HCC patients, we assumed that CAGS may be helpful in predicting the response to ICB immunotherapy. After categorizing the cancer patients based on the optimal cutoff value of CAGS, those with a high CAGS exhibited a clinical benefit in both the anti‐PD‐1 (GSE78220, Fig. 3A,B) and anti‐PD‐L1 (IMvigor210, Fig. 3C,D) ICB therapy cohorts. Furthermore, as seen in Fig. 3E, the high CAGS group had a higher responsive rate (response/nonresponse = 48.33%/51.67%) to ICB treatment than the low CAGS group (response/nonresponse = 16.39%/83.61%). However, with the combination of CAGS and tumor neoantigen (NEO) burden, it was observed that the CAGS‐high/NEO‐high group had a comparable survival probability as the CAGS‐low/NEO‐high group, which suggested that the survival impact of CAGS may be inferior to that of NEO level (Fig. 3F).

Fig. 3

Prediction of immune checkpoint blockage therapy response. To detect the potential ability of CAGS in predicting ICB treatment response, the CAGS were quantified using IMvigor210 urothelial cancer and GSE78220 melanoma ICB cohorts. Patients with high CAGS exhibited a better overall survival condition than that with low CAGS in both the anti‐PD‐1 (GSE78220, A–B) and anti‐PD‐L1 (IMvigor210, C–D) ICB therapy cohorts. (E) The high CAGS group had a higher responsive rate (response/nonresponse = 48.33%/51.67%) to ICB treatment than the low CAGS group (response/nonresponse = 16.39%/83.61%). (F) Overall survival analysis of ICB therapy cohorts based on both CAGS and NEO status. CAGS, cohesin‐associated gene scores; ICB, immune checkpoint blockage; NEO, neoantigen burden.

Comprehensive overexpression of STAG1 in global HCC mRNA data and in‐house IHC validation

Because STAG1 is the core subunit of the cohesin complex protein and it serves as a TF, its expression status in HCC tissue samples was further explored. A total of 37 platform matrices were enrolled in this study (Table S1), which covered 76 datasets from 12 countries (Table S2). The SMD forest plot indicated that STAG1 was upregulated in 3313 HCC tissue samples in comparison with 2692 non‐HCC tissue samples (Fig. S4A). The sensitivity analysis plot indicated that the included datasets could not explain the major source of heterogeneity (Fig. S4B). The funnel plot implied insignificant publication bias, which showed the stability of the quantitative synthesis result (Begg's test: P value = 0.067) (Fig. S4C). The summary characteristics operating curve showed a moderate discriminatory ability of STAG1 (Fig. S4D), with weak sensitivity (Fig. S4E) and moderate specificity (Fig. S4F). Fagan's nomogram and likelihood ratio forest plots indicated the general accuracy of STAG1 in differentiating between the HCC and non‐HCC tissue samples (Fig. S4G–I). However, the OS prediction ability of STAG1 was insignificant in HCC patients (data not shown).

More importantly, the STAG1 protein was remarkably stained in the nucleus of the HCC cells (Fig. 4A–D) in comparison with the non‐HCC cells (Fig. 4E,F), thus validating the overexpression trend of STAG1 in HCC tissue samples (Fig. 4G,H).

Fig. 4

In‐house immunohistochemistry staining of STAG1 transcriptional factor in HCC tissue samples. STAG1 protein was stained in the nucleus of the HCC cells (A–D) and compared with healthy cells (E–F). Increased protein expression of STAG1 in HCC displayed a strong ability in differentiating HCC from non‐HCC tissue samples (G–H). HCC, hepatocellular carcinoma. The scale bars in panels A–F represent 100, 50, and 20 μm, from left to right.

Prospective signaling pathways and transcriptional regulation of STAG1 in HCC

By adopting a network embedding algorithm, a total of five predominant functional modules were identified from the downregulated STAG1 negative CEGs (Fig. 5), which were enriched in lymphocyte activation (c1_3), monocarboxylic acid metabolic process (c1_4), lipid catabolic process (c1_7), cell junction organization (c1_8), and blood vessel development (c1_38). Additionally, we constructed a scale‐free co‐expression network based on the upregulated STAG1‐positive CEGs (soft threshold = 19) (Fig. 6A). The blue module, defined as the high‐CAGS HCC phenotype‐related module (Fig. 6B,C), was predominantly enriched in the mitotic cell cycle process (Fig. 6D). The similarity scores of the STAG1 putative transcriptional targets were calculated, and PDS5 cohesin‐associated factor A (PDS5A) and inhibitor of DNA binding 1, HLH protein (ID1) were identified as the core genes in the mitotic cell cycle process and blood vessel development, respectively (Fig. 6E). More importantly, the downregulated platelet‐derived growth factor receptor alpha (PDGFRA) in blood vessel development (Fig. 7A) and upregulated PDS5A in the mitotic cell cycle process (Fig. 7B) were predicted as two transcriptional targets of STAG1.

Fig. 5

Potential functional modules of downregulated STAG1 negative co‐expression genes in HCC. Based on the downregulated STAG1 negative co‐expressed genes, a total of five predominant modules were identified using the multiscale embedded gene co‐expression network analysis (MEGENA).

Fig. 6

Prospective signaling pathways of upregulated STAG1 positive co‐expression genes in hepatocellular carcinoma. We performed weighted gene co‐expression network analysis (WGCNA) to determine the high‐CAGS HCC phenotype‐related module genes. (A) A soft threshold of 19 was set to construct a scale‐free topological model. (B–C). A blue gene module exhibited a highly positive correlation with the CAGS‐high phenotype in HCC. (D) The blue gene module was predominantly enriched in the mitotic cell cycle process. (E) The similarity scores of the STAG1 putative transcriptional targets were calculated, and PDS5A and ID1 were identified as the core genes in blood vessel development and mitotic cell cycle process, respectively. CAGS, cohesin‐associated genes scores; HCC, hepatocellular carcinoma.

Fig. 7

Transcriptional regulation of STAG1 in mitotic cell cycle process and blood vessel development. We identified two transcriptional targets of STAG1 in hepatocellular carcinoma. (A) PDS5A in the mitotic cell cycle process. (B) PDGFRA in blood vessel development.

As is shown in Fig. S5A, PDGFRA mRNA was significantly downregulated in 3348 HCC when compared with 2954 normal liver tissue specimens. Decreased PDGFRA showed a moderate discriminatory ability between HCC and normal liver tissue (Fig. S5B,C). The STAG1 expression level was inversely correlated with the mRNA level of PDGFRA (Fig. S6). Based on the potential role that STAG1 has on blood vessel development in HCC, we subsequently explored whether STAG1 could transcriptionally regulate the transcription of vascular endothelial growth factor A (VEGFA). Intriguingly, STAG1 showed an obvious transcriptional factor binding intensity in the promoter region of VEGFA (Fig. S7A). Furthermore, STAG1 was positively correlated with the expression level of VEGFA (Fig. S7B).

Promising anti‐HCC agents by targeting STAG1 transcriptional mechanisms

Based on the enriched gene ontology terms of blood vessel development and mitotic cell cycle process, a total of 91 putative transcriptional targets were inputted to screen for sensitive anti‐HCC agents, either FDA‐approved drugs or clinical trial drugs, and were used to identify the potential interventional targets for treating HCC patients. Surprisingly, downregulated endothelial PAS domain protein 1 (EPAS1) (SMD = −0.6386, P value < 0.0001) exhibited a significant negative correlation with the IC50 values of 14 agents (such as TAK‐960 analog, methotrexate, benzimate, TAK Plk inhibitor, oxaliplatin, volasertib, and teglarinad) and was positively correlated with that of XAV‐939 and AZD‐8055 (Figs S8 and S9). Among them, daporinad (a nicotinamide phosphoribosyl transferase inhibitor), EMD‐534085 (a mitotic kinesin‐5 inhibitor), GSK‐461364 (a polo‐like kinase 1 inhibitor), and teglarinad (a prodrug of nicotinamide phosphoribosyl transferase inhibitor) ligands could be docked by the EPAS1 protein with a high total score (all with total scores > 5.0). Further interaction analysis implied that the hydrogen bond was the predominant type of force (Fig. 8).

Fig. 8

Putative therapeutic agents for treating hepatocellular carcinoma. Molecular docking between EPAS1 and four antitumor agents, (A) Daporinad (nicotinamide phosphoribosyl transferase inhibitor). (B) EMD_534085 (mitotic kinesin‐5 inhibitor). (C) GSK_461364 (polo‐like kinase 1 inhibitor). (D) Teglarinad (a prodrug of nicotinamide phosphoribosyl transferase inhibitor), which has been approved by the US Food and Drug Administration.

Discussion

HCC is a lethal abdominal malignancy with few effective intervention strategies in its late stages [37, 38]. The advent of immunotherapy, represented by ICB therapy, created novel opportunities for advanced HCC patients [39, 40, 41]. In this study, we established a CAGS‐based signature scoring system and identified three CAG patterns in HCC patients, which may be used to define distinct clinical prognosis, TME, and genetic mutation characterization. Moreover, a high CAGS may be a potential indicator for therapeutic response to anti‐PD‐1/PD‐L1 ICB treatment. STAG1, as an important CAG TF, was comprehensively explored in the transcriptional regulatory mechanisms of blood vessel development and the mitotic cell cycle process in HCC, where EPAS1 may be a drug‐treatment target for HCC treatment.

To date, numerous biomarkers and models have been found to predict immunotherapy response [42, 43, 44]. The novelty of the present study is that we predicted the response status to ICB therapy by converting 17 CAG signatures into personalized ssGSEA scores [45], rather than by directly depending on the expression values. Intriguingly, a high CAGS predicted poor OS prognosis in HCC patients but exhibited a clinical advantage in anti‐PD‐1/PD‐L1 immunotherapy, with a higher response rate in comparison with the low CAGS group (response/nonresponse = 48.33%/51.67%). Additionally, a combination of high NEO burden [46] and high CAGS may be used to predict a better clinical prognosis when receiving standard ICB treatment.

It is well‐known that TME plays a pivotal role in promoting the invasion and metastasis of HCC by regulating immune and stromal cells [47, 48]. In this study, we identified three CAGS patterns from the HCC patients and revealed the distinct clinical and molecular characteristics among them. We noticed that the CAGS cluster 2 group displayed an attenuated toxic metabolism level and had an increased mutation rate of TP53. Further exploration revealed that the CAGS cluster 2 was characterized by high tumor purity along with deceased immune and stromal scores. To the best of our knowledge, cohesin mutations are commonly detected in a range of malignancies [49, 50, 51], and the TP53 suppressor gene is one of the most common mutation genes in HCC [52]. More importantly, there are intimate associations between cohesin and TP53 mutations. For instance, there is a high co‐occurrence between TP53 gene alterations and cohesin mutations [53]. The STAG1 cohesin subunit has been reported to serve as a direct transcriptional target of TP53 [54]. In response to genotoxic stresses, SATG1 mRNA expression could be induced in a TP53‐dependent manner, and endogenous STAG1 protein expression increased significantly when exposed to ultraviolet radiation [54]. Notably, STAG1 also participated in DNA damage repair [55], in addition to the renowned genome guardian, TP53 [56]. Moreover, TP53 gene mutation abrogated its role in promoting DNA damage repair and cellular apoptosis, thus leading to severe cell carcinogenesis [57]. In this setting, CAGs may be used to define a poorer HCC phenotype with higher tumor purity and more TP53 mutations.

Furthermore, the present study shed light on the transcriptional mechanisms of the STAG1 cohesin subunit involving blood vessel development and the mitotic cell cycle process in HCC. We first confirmed the overexpression of STAG1 in 3313 HCC tissue samples based on multicentered HCC bulk RNA‐seq data (SMD = 0.54). Intriguingly, it was observed that the expression level of STAG1 was increased after hsa‐miR‐23a or hsa‐miR‐27a knockout [58], which suggested that microRNA was involved in the overexpression of STAG1 in HCC. Furthermore, according to six ChIP‐seq datasets [59, 60, 61], the PDS5A and PDGFRA promoter regions both exhibited remarkably high TF‐DNA binding intensities of STAG1 (range: 0–20). Based on the global HCC mRNA expression data from 12 ethnicities, PDS5A was an upregulated gene and positively correlated with STAG1; conversely, PDGFRA was downregulated and negatively correlated with STAG1.

Similar to the cohesin core subunit STAG1, the PDS5A regulatory subunit is mainly involved in modulating sister chromatid cohesion during mitosis [62]. The PDS5A protein participates in the dynamic association between cohesin and chromatin and protects the replication forks from degradation [63]. After the deletion of the PDS5 protein, cell cycle progression was inhibited [64]. In this context, STAG1 TF is proposed to positively regulate cell cycling by targeting PDS5A, thus promoting mitotic division and malignant proliferation of HCC cells. Additionally, STAG1 was correlated with vascular development. It was demonstrated that STAG1 and STAG2 balanced the production of hematopoietic/vascular progenitors in a Zebrafish embryo model [65]. In the shTMEM30A‐treated primary human retinal endothelial cells, STAG1 was significantly downregulated and the vascular formation was inhibited [66]. Herein, we proposed that PDGFRA may be a transcriptional target of STAG1 in HCC. PDGFRA encodes a tyrosine kinase receptor and is involved in embryonic development and wound healing [67]. The PDGFRA signaling pathway may promote angiogenesis, carcinogenesis, and tumor dissemination [68, 69]. Moreover, PDGFRA was downregulated in the HCC tissue samples but overexpressed in endothelial cells, and it correlated with vascular invasion and forecast poorer OS and disease‐free survival in HCC patients [70]. Taken together, these studies showed that downregulated PDGFRA may be negatively regulated by STAG1 TF in HCC cells but positively correlated with the emerging vessels surrounding the tumor lesion. To support our hypothesis, we also investigated the transcriptional activity of STAG1 and the tumor angiogenesis promoting factor, VEGFA. Surprisingly, STAG1 was positively correlated with VEGFA and exhibited a specific transcriptional factor binding peak at the upstream promoter region of the VEGFA gene. It was conceivable that STAG1 may promote the angiogenesis in HCC by positively regulating the transcription of VEGFA. The present study revealed two important transcriptional targets for STAG1 in HCC development. In the future, in vitro and in vivo assays must be conducted to validate the transcriptional mechanism of STAG1 underlying HCC.

Several anti‐HCC agents were also identified by targeting the EPAS1 protein, all with high total scores. Interestingly, a previous study had reported that microRNA‐3609 could attenuate the expression of EPAS1 and sensitize HCC cells to sorafenib treatment [71]; that study was the first to reveal the potential role of EPAS1 in HCC treatment. The present study further expanded on that finding by identifying EPAS1 as the targetable protein in treating HCC. Specifically, teglarinad (GMX1777) exhibited the highest total score when docked by the EPAS1 protein, and it was found that higher EPAS1 expression correlated with higher sensitivity of teglarinad to treat HCC cells. Teglarinad is believed to suppress the development of tumor cells by interfering with DNA repair, abrogating angiogenesis, and enhancing radiation efficacy [72]. However, it is necessary to verify the potential implications of the screened anti‐HCC agents at the cell and animal levels, or even at the clinical trial level. It is believed that the development of emerging technologies of in‐silico drug designs and artificial intelligence will result in more effective therapies for curing HCC patients in the future.

Conclusions

In this study, we identified three distinct CAG phenotypes from HCC patients, and a high CAGS predicted poor prognosis of HCC patients but correlated with potential clinical advantages in ICB immunotherapy. STAG1 may fuel HCC deterioration by transcriptionally regulating blood vessel development and the mitotic cell cycle.

Conflict of interest

The authors declare no conflict of interest.

Author contributions

(I) Study design, experimental supervision, and manuscript editing: C‐ZL, GC, R‐QH, Z‐GH, J‐JL, X‐FD, and X‐PL; (II) Experimental operation, data analysis and interpretation, and manuscript writing: C‐ZL, J‐DL, and RL; (III) Final approval of manuscript: All authors.

Fig. S1. Design route of this research.

Fig. S2. Metabolic characterization and mutation co‐occurrence statuses of HCC patients. Based on the identified CAG clusters, the metabolic characterization and mutation co‐occurrence statuses of different HCC patients were compared. (A–C) The HCC patients in CAG cluster 2 with poor prognosis displayed prominently lower metabolism levels of drug, retinol, and xenobiotics than the HCC patients in the other CAG clusters. Panels D–F show the mutation co‐occurrence statuses of different CAG patterns. TP53 mutation co‐occurrence was more frequent in the HCC patients in the CAG cluster 2 than in the patients in the other CAG clusters. Abbreviations: HCC, hepatocellular carcinoma; CAG, cohesin‐associated gene.

Fig. S3. Prognostic prediction ability of STAG1 in patient receiving immunotherapy. We appraised the prognostic value of STAG1 in the IMvigor210 urothelial cancer (A) and GSE78220 melanoma (B) immunotherapy cohorts. Higher STAG1 expression level predicted better overall survival outcomes in patients who received immunotherapy.

Fig. S4. Upregulation of STAG1 based on global HCC data. A total of 37 platform matrices were enrolled to analyze the overall expression status of STAG1 in HCC. (A) The standard mean difference forest plot indicated that STAG1 was upregulated in 3313 HCC tissue samples in comparison to 2692 non‐HCC tissue samples. (B) The sensitivity analysis plot indicated that the included datasets could not explain the major source of heterogeneity. (C) The funnel plot implied insignificant publication bias, which showed the stability of the quantitative synthesis result (Begg's test: P value = 0.067). (D) The summary characteristics operating curve showed a moderate discriminatory ability of STAG1, with weak sensitivity (E) and moderate specificity (F). (G–I) Fagan's nomogram and likelihood ratio forest plots indicated the general accuracy of STAG1 in differentiating between HCC and non‐HCC tissue samples. Abbreviation: HCC, hepatocellular carcinoma.

Fig. S5. Comprehensive mRNA expression level of PDGFRA in the HCC tissue samples. PDGFRA was significantly downregulated in the HCC tissue samples when compared with normal liver tissue specimens (A). Downregulated PDGFRA mRNA showed a moderate discriminatory ability between HCC and normal liver tissue samples (B, C). Abbreviation: HCC, hepatocellular carcinoma.

Fig. S6. Negative correlations between STAG1 factor and PDGFRA target. We computed the Pearson correlation coefficients between STAG1 expression level and PDGFRA expression level. STAG1 expression level was inversely correlated to PDGFRA mRNA level.

Fig. S7. Potential correlations between STAG1 factor and the VEGFA mRNA expression. Chromatin immunoprecipitation sequencing data were reanalyzed to explore the potential regulation between STAG1 and vascular endothelial growth factor A (VEGFA). STAG1 showed an obvious transcriptional factor binding intensity in the promoter region of VEGFA (A). STAG1 was positively correlated to the expression level of VEGFA (B).

Fig. S8. Potential correlations between EPAS1 expression and the sensitivity of anti‐cancer agents. We computed Pearson correlation coefficients between EPAS1 expression and half maximal inhibitory concentration of the screened agents.

Fig. S9. The discrepancies of half maximal inhibitory concentration in the EPAS1‐high expression group and the EPAS1‐low expression group. Note: A Wilcoxon test was conducted to compare the discrepancies of half maximal inhibitory concentration between EPAS1‐high expression group (n = 30) and the EPAS1‐low expression group (n = 30). ***P value <0.001.

Table S1. Overexpressed STAG1 in 37 hepatocellular carcinoma platform matrices.

Table S2. Fundamental dataset information of the included global hepatocellular carcinoma datasets.

Click here for additional data file.