Development and Validation of a Novel Immune-Gene Pairs Prognostic Model Associated with CTNNB1 Alteration in Hepatocellular Carcinoma

Background

Despite the great progress made in treatment of hepatocellular carcinoma (HCC), the prognosis for patients with the disease remains poor [1–3]. In recent years, immunotherapies and particularly immune checkpoint inhibitors, has revolutionized treatment of many malignancies. HCC is a typical inflammation-related cancer that occurs in the presence of immunosuppression, making immunotherapy a potentially attractive treatment option [1,4–6]. However, some recent clinical trials have found that patients with CTNNB1 mutations in HCC tissues are not sensitive to immunotherapy [7–9]. The underlying mechanism is not entirely clear.

CTNNB1 is an important component of the typical Wnt signaling pathway and plays an indispensable role in physiological embryology, zone formation and metabolic regulation of the liver [10]. Activating and inactivating mutations in linker b1 (CTNNB1), which encodes the protein β-catenin, are associated with tumorigenesis of the liver in humans and [11] Approximately 11% to 41% of all liver malignancies are HCC with mutated CTNNB1,which has specific metabolic, morphological, and clinicopathological characteristics [12]. For example, tumors with CTNNB1 mutations often occur in the setting of a liver without cirrhosis, the infection rate of hepatitis b virus is low, and the differentiation of the tumor is good, but it is often accompanied by microvascular invasion, tumor capsule invasion, and intratumor cholestasis [13]. Mutations in CTNNB1 also are associated with malignant transformation of hepatocellular adenoma [10].

A recent study found that activation of β-catenin promoted immune escape as well as resistance to anti-programmed cell death protein-1 therapy for HCC. The mechanism may be related to defective dendritic cells (DC) and antigen-specific T-cell recruitment [9]. In a mouse model, activation of β-catenin led to poor T-cell infiltration, promoting immune escape and leading to resistance to immunotherapy.7 Although the relationship between CTNNB1 mutations and prognosis remains controversial, there is increasing evidence suggesting that changes in the tumor immune response can affect prognosis in patients with HCC [14,15].

In this study, we attempted to explore the immunophenotype of HCC from the perspective of mutations in CTNNB1. We identified differentially expressed immune genes (DEIGs) in samples from patients with and without CTNNB1 mutations in The Cancer Genome Atlas (TCGA) database and paired them to explore any correlation with prognosis. Creation of our prognostic model was based on the above DEIGs. In addition, the International Cancer Genome Consortium (ICGC) dataset was used to validate the accuracy and universal applicability exhibited by the prognostic model.

Material and Methods

DATA COLLECTION:

mRNA Sequence data from 374 HCC samples came from TCGA (https://portal.gdc.cancer.gov/repository). Corresponding clinical data – including survival time, survival status, sex, age, race, body mass index (BMI), histological grading, alpha fetoprotein (AFP) level, vascular invasion, American Joint Committee on Cancer tumor, node, metastasis (AJCC-TNM) stage, family history of cancer, prior malignancy history, and new tumors presenting following initial treatment, as well as tumor status – were taken from the University of California Santa Cruz’s Xena Funcational Genomics Exporer (https://xenabrowser.net/). The CTNNB1 mutation sample list was obtained from the cBioPortal (https://www.cbioportal.org/). We downloaded another independent mRNA sequence dataset and corresponding clinical data from the ICGC (https://icgc.org/) to conduct external validation. The two sequence datasets came from the Illumina HiSeq RNA-Seq platform. We downloaded the immune-related gene list from the ImmPort database (https://immport.niaid.nih.gov). It should be stated that acquisition of the above data followed the data access policies and release guidelines for these databases. Therefore, there was no requirement to obtain the approval from local ethics committee.

IDENTIFICATION OF DEIRGS:

The “edgeR” R package was used to identifying differentially expressed immune-related genes (DEIRGs) in 276 HCC samples free of CTNNB1 mutations and 96 HCC samples with CTNNB1 mutations in the TCGA cohort. A false discovery rate (FDR) of <.05 was considered significant.

ANNOTATIONS OF DEIRG BIOLOGICAL PROCESSES:

The “clusterProfiler” R package was used to annotate the function of the DEIRGs and conduct enrichment analysis of their involvement in biological processes (gene ontology). The threshold value was an FDR <.05.

CONSTRUCTION AND VALIDATION OF IMMUNE GENE-PAIRS PROGNOSTIC MODEL:

We then paired the DEIRGs. In each gene pair, if the expression level of the former gene was higher than that of the latter gene, the value assigned was 1. On the contrary, if the expression level of the former gene was lower than that of the latter gene, the value assigned was 0. Immune-gene pairs with values of 0 or 1 that comprised less than 20% of the TCGA cohort were excluded. A total of 14 485 immune gene pairs were obtained after screening. Univariate Cox regression analysis was applied to identify the immune gene pairs associated with prognosis (P<.001). The “|glmnet” R package was used for Least absolute shrinkage and selection operator (Lasso) regression analysis of the prognostic immune gene pairs. The penalty coefficient (log Lambda value) with the smallest cross-validation error helped to confirm the gene pair coefficient. The formula for the risk score was the sum of each coefficient of immune-gene-pair multiplied by each value in the immune-gene pair. The risk score corresponding to the maximum value for sensitivity and specificity under the receiver operating characteristic (ROC) curve was used to classify the groups with high and low risks. The R packages “survminer” and “survivalROC” were used to generate Kaplan-Meier survival curves and ROC curves regarding the risk score as a means of assessing the predictive power of the model. The log-rank test was used to compare the two groups with regard to survival curves, and P<.05 was considered statistically significant. Patients in TCGA were divided into several subgroups based on their clinical characteristics for internal verification, and external verification was conducted by using the ICGC cohort.

INDEPENDENCE VALIDATION OF THE PROGNOSTIC MODEL:

Univariate and multivariate Cox regression analysis assisted in verifying whehter risk score could independent predict HCC prognosis in the TCGA and ICGC cohorts.

ANALYSIS OF ASSOCIATION BETWEEN RISK SCORE AND CLINICOPATHOLOGY:

We used the Wilcoxon test or Kruskal’s algorithm (multiple groups) to analyzing the association between the risk score and clinicopathology (including AJCC-TNM stage, histologic grade, vascular tumor cell type, tumor status, and new tumor event after initiation of treatment). P<0.05 was deemed to be statistically significant. The “beeswarm” package in the R software was used to generate the boxplot.

ANALYSIS OF INFILTRATION OF IMMUNE CELLS IN DIFFERENT RISK GROUPS:

The CIBERSORT algorithm assesses cell composition in complex tissues based on standardized gene presentation data, quantifying the richness of specific kinds of cells [16]. We used this method (with a characteristics matrix of 547 genes standing for 22 kinds of infiltrating immune cells) to quantify immune cell infiltration in tumor tissues from groups with high and low risks, respectively. P<.05 was used as the standard by which to judge accuracy of the estimation.

STATISTICAL ANALYSIS:

Statistical analyses were conducted with R software v3.6.1 (R Foundation for Statistical Computing, Vienna, Austria) as well as with GraphPad Prism v7.00 (GraphPad Software Inc., United States). Analysis of qualitative variables was performed with Pearson’s χ² test or Fisher’s exact test; analysis of quantitative variables was done with a non-parametric Wilcoxon test for unpaired samples in the appropriate mode. Analysis of various normalized data groups was done with Kruskal’s algorithm. P<.05 was considered statistically significant.

Results

IDENTIFICATION OF DIFFERENTIALLY EXPRESSED IMMUNE GENES:

Immune cell infiltration in tumor tissues in the wild-type (WT) CTNNB1 group was significantly different from that in the mutated CTNNB1 group (Figure 1A). Among the 434 differentially expressed immune genes (the heat map shown in Figure 1B), 231 were upregulated in the non-CTNNB1-mutated group (n=276), and their biological processes mainly involved activation and regulation of leukocytes and lymphocytes (Figure 1C). We found that 203 genes were upregulated in the CTNNB1-mutated group (n=98), and their biological processes were mainly related to innate immune response and to regulation of the Fc receptor signaling pathway. (Figure 1D). This may help to explain the insensitivity of CTNNB1-mutated HCC to immunotherapy.

CONSTRUCTION OF A ROBUST PROGNOSTIC MODEL CONSISTING OF 45 IMMUNE GENE PAIRS:

Univariate Cox regression analysis showed that among the 14 485 gene pairs, 480 i were associated with prognosis (P<.001). After the overfitting was removed using Lasso regression, a robust prognostic model consisting of 45 immune gene pairs was established (Table 1). According to the ROC curve, the optimal critical value for dividing the group in to highand low risk was.886. Kaplan-Meier survival curve analysis showed that the group with high risk had lower rates of overall (OS) and disease-specific survival (DFS) relative to the group with low risk (Figure 2A, 2C). AUC values for 1-, 2-, 3-, 4-, and 5-year OS reached.899,.905,.892,.915, and.894, respectively (Figure 2B) and those for DFS were.956,.945,.936,.935, and.903, respectively (Figure 2D). As displayed in the survival status map, a risk score increase meant an obvious increase in the number of patient deaths (Figure 3A). As illustrated in univariate and multivariate Cox regression analyses, only risk score was an independent predictor of prognosis (Figure 3B, 3C).

INTERNAL VALIDATION BASED ON CLINICAL FEATURES OF TCGA:

To observe whether the model constructed was applicable to different populations, we divided the patients into 29 subgroups according to their clinical characteristics (including sex, age, race, BMI, alpha-fetoprotein levels, histological grade, vascular invasion, AJCC-TNM stage, family history of cancer, malignancy history, new tumors presenting following initial treatment, and individual tumor status) for internal verification. The survival curve showed that in each subgroup, the group with a high risk had lower OS relative to the group with low risk (Figure 4A–4L). This indicated that our prognostic model had universal applicability.

CORRELATION ANALYSIS BETWEEN RISK SCORE AND CLINICOPATHOLOGY:

With increasing histological grade, AJCC stage, and vascular invasion, patients’ median risk scores gradually increased. Median risks score in patients with new tumors after initial treatment were higher than in those without new tumors, and median risk scores in patients with a tumor was higher than those in individuals without tumors (Figure 5). We also carried out Chi-squared tests on the two groups, which revealed an obvious difference in clinical stage, histological grade, vascular invasion degree, and tumor histological grade between them, which may help to explain the differences in their prognoses (Tables 2, 3).

EXTERNAL VALIDATION REGARDING THE PROGNOSTIC MODEL IN THE ICGC COHORT:

Patient risk scores in the ICGC cohort were calculated based on the risk score calculation formula from the TCGA cohort. The cutoff values of the groups with high and low risks were the same as those for the TCGA cohort. As with the TCGA cohort, the group with a high risk had a poor prognosis (Figure 6A). AUC values for 1-, 2-, 3-, 4-, and 5-year OS reached.831,.727,.730,.786 and.786, respectively (Figure 6B). As revealed by univariate and multivariate Cox regression analyses, risk score was an independent predictor of prognosis (Figure 6C, 6D).

INFILTRATION OF IMMUNE CELLS IN TUMOR TISSUE BETWEEN DIFFERENT RISK GROUPS:

After screening (P<.05), 42 tumor tissues from the low-risk group and 38 tumor tissues from the high-risk group were included. The infiltration ratio for B cells, CD8 T lymphocytes, activated natural killer cells, and M1 macrophages in the group with low risk appeared higher relative to the group with high risk, while the latter group had a higher infiltration ratio for M0 and M2 macrophages relative to the former group (Figure 7).

Discussion

The prognosis for patients with HCC is still not satisfactory, mainly because most cases are only diagnosed when the disease is advanced [1]. Adjuvant therapy, such as immunotherapy, is one of the primary treatments for advanced HCC. However, several recent clinical trials have found that patients who have HCC with CTNNB1 mutations are not sensitive to immunotherapy [5,8,9]. The mechanism of CTNNB1 mutation affecting immunophenotypic regulation and prognosis of HCC remains unclear. Some immune signatures related to prognosis of HCC have been identified. For example, Long [17] identified a dual-immune-gene prognostic model, Wang [18] identified a 9-immune-gene prognostic signature, and Li [19] identified a 6-immune-gene prognostic signature. These reports provide theoretical evidence on which to base assessment of HCC prognosis, taking into account genes related to the immune system.

In this study, we found that several immune cells in tumor tissues with CTNNB1 mutation had a significantly decreased infiltration level, which is similar to previous reports [9,15]. In HCC without CTNNB1 mutations, the upregulated immune genes were significantly involved in regulation of lymphocyte activation. Therefore, immunotherapy resistance associated with a CTNNB1 mutation may be related to low lymphocyte activity.

After it was clear that the differentially expressed genes could significantly affect the immune response of the tumor, we combined them in pairs. This is a different approach than that in previous studies. Most previous studies have assessed the association between a single gene and prognosis based on its level of expression. With gene pairing, we were able to analyze the effect of the relative expression level of two genes on prognosis, transforming continuous variables into a dichotomous variable. With this methods, other centers using the model would not have to carry out batch normalization.

We screened 45 gene pairs from 480 immune gene pairs related to prognosis by means of lasso regression. The results showed that risk stratification of patients with based on these 45 immune gene pairs was very accurate, especially for evaluating the specific survival rate for HCC. It was important that the prognostic model had universal applicability and could be applied to people with different clinical characteristics. In addition, risk scores were associated with factors correlated with poor prognosis, such as low tumor differentiation, vascular invasion, and advanced clinical stage. This may help to explain the difference in prognosis between different groups.

Finally, we used the ICGC dataset to conduct external verification of the prognostic model. Conforming to the findings of the TCGA cohort, the prognostic model was highly accurate for predicting survival in the short and long term, and it was an independent predictor of prognosis. We compared AUC values to determine the accuracy of the prognostic model relative to other prognostic signatures. AUC values for Long’s double immune gene prognostic signature were.738,.727, and.655 at 1, 3, and 5 years, respectively [17]. AUC values for Wang’s 9-immune-gene prognostic signature were.811,.711, and.734 at 1, 3, and 5 years, respectively [18]. AUC values for Li’s 6-immune-gene prognostic signature were.761,.681, and.692 at 1, 3, and 5 years, respectively [19]. AUC values for our prognostic model at 1, 3, and 5 years were.899,.892, and.894, respectively. Obviously, our prognostic model has more value than any of the traditional prognostic immune signatures.

Tumor microenvironment (TME) is the cellular environment of tumor cells, and includes tumor stromal cells, the extracellular matrix, and soluble molecules. Once the TME is shaped, a lot of immune cells, such as T cells, macrophages, and myelogenic inhibitory cells, are drawn by chemotaxis to constitute it [20]. In the TME, immune and stromal cells, which are major types of non-tumor components, may have a role in diagnosis and prognostication of tumors [21,22]. As found in this study, prognosis for patients in the group with low risk was dramatically better than in the group with high risk. Garnelo [23] found that patients with higher proportions of tumor-infiltrating T and B cells had stronger local immune activity and better prognosis. This was consistent with our research results. In addition, patients whose activated NK-cell infiltration level was high had better prognosis. That is in accord with conclusions by Zhu and other researchers [24].

Tumor-associated macrophages (TAMs) are macrophages that infiltrate tumor tissues and they make up the largest number of immune cells in the TME [25,26]. M1 macrophages have an antitumor effect, which can distinguish tumor cells from normal cells and eventually kill tumor cells by identifying them. M2 macrophages mainly promote tumor development, invasion, and metastasis [27,28]. In patients who had a high level of infiltration of M2 macrophages, the prognosis was poor, while for those with high levels of M1 macrophage infiltration, the prognosis was better. Therefore, regulating the phenotype of TAM may function as a positive factor in improving prognosis in patients with HCC.

Some genes in our model have been proven to be related to tumor prognosis and immune regulation. For example, AHNAK can promote tumor metastasis via transforming growth factor-β-mediated epithelial-mesenchymal transition [29], and CCL21 can promote immune activity in the TME through co-localization of dendritic cells (DC) and play a role in T-cell programming of ectopic lymph node architectural structures, which has been correlated to cancer prognosis [30]. CD48 is indispensable as an environmental sensor for regulating HSC and progenitor cell numbers, and it also inhibits tumor development [31]. Differential activation of the transcription factor IRF1 laid the foundation for the distinct immune responses elicited by Type I and Type III interferons [32]. Levels of Januse kinase 1 mRNA have been correlated with prognosis and degree of immune infiltration of breast cancer [33]. Matrix metalloproteinase-12 (MMP-12) over-presentation is associated with poor prognosis in HCC [34]. Stimulating MSR1 can enhance c-Jun N-terminal kinase -mediated inflammation in interleukin-4-activated macrophages [35]. NEDD4 exerts a key effect in enhancing HCC multiplication and metastasis by activating the PTEN/phosphatidylinositol-3 kinase/Akt signaling pathway [36]. SEMA3F can accelerate HCC metastasis by activating the focal adhesion pathway [37].CCL5 is a mediator for CD40-driven CD4+ T-cell tumor infiltration as well as immunity in pancreatic cancer.

Conclusions

Our study provides a new way to evaluate prognosis in patients with HCC, but there are still some limitations. Although our prognostic model has been verified internally and externally, this was a retrospective study, so further prospective clinical trials are necessary. The relative expression level exhibited by the two immune genes had a significant impact on prognosis and further study of its mechanism is warranted.

This research suggests a novel approach for evaluating the prognosis of HCC, which may be beneficial for guiding individualized treatment of patients.