Construction and validation of a novel IGFBP3-related signature to predict prognosis and therapeutic decision making for Hepatocellular Carcinoma

Public data acquisition and processing

The transcriptome data and relevant prognostic resource of HCC in the Cancer Genome Atlas (TCGA, https://portal.gdc.cancer.gov/repository), International Cancer Genome Consortium (ICGC) Japanese liver cancer (ICGC-LIRI-JP) cohort (https://dcc.icgc.org/projects/LIRI-JP), GSE54236, GSE14520 (GPL3921), and GSE76427 were downloaded and processed as reported publications (Zhang et al., 2022a; Zhang et al., 2022b). For analysis, we applied log₂[TPM+1] transformed expression data. IGFBP3 protein levels of differentiation analysis with the aid of CAPTC database (https://cprosite.ccr.cancer.gov/). In brief, select “Liver Cancer” from the drop-down menu of “Tumor Types”, “IGFBP3” from the drop-down menu of “Gene” on the page, click “Submit” for analysis, and then click “Export Data” to obtain protein expression profile data. Figure 1 shows the analysis process of this study. The data sources and related sample numbers are detailed in Table S1.

Cell culture and RT-qPCR

Human normal liver cell line (LO2) and human liver cancer cell line (Huh-7) were purchased from Hongshun Biotechnology Co. LTD (Shanghai, China). LO2 cells were cultured in RPMI-1640 (Procell, Wuhan, China) containing 20% FBS (Procell, Wuhan, China). Huh-7 cells were cultured in DMEM (Gibco, Billings, MO, USA) with 10%FBS. All cells were cultured in a 5% CO2 incubator humidified at 37 °C. Total RNA was isolated using SteadyPure Quick RNA Extraction Kit (AG21023; Accurate Biology, Changsha, China) according to the manufacturer’s manual. cDNA was synthesized by MCE RT Master Mix for qPCR II (MCEs, NJ, USA). A GoTaq® qPCR Master Mix (A6001; Promega) was used for qPCR. Primers were synthesized by Shangya Biotechnology (Fuzhou, China). IGFBP3: forward: 5′-AAATGCTAGTGAGTCGGAGGA- 3′, reverse: 5′-CTCTACGGCAGGGACCAT ATT- 3′. We used GAPDH as a reference for IGFBP3 following the 2^−ΔΔCT method.

Tumor tissues specimens and immunohistochemistry

Tissue specimens were fixed by 10% formalin and embedded in paraffin. The tissues were then cut into 3um thick sections. After dewaxing and hydration, citric acid buffer (0.01M, pH 6.0) was used and boiled for 15 min for antigen repair. Immunohistochemically staining was performed using the EliVision™Plus kit (Maixin Biotechnology, Fuzhou, China). Subsequently, sections were incubated overnight with anti-IGFBP3 polyclonal antibody (ER1911-12) (1:200) or PBS (negative control) at 4 °C, then coupled with secondary antibody at room temperature for 10 min and stained with diaminobenzidine (DAB Kit, Lab Vision) for 40 s. The cells of the patients with positive immunohistochemically reaction were stained with hematoxylin for 15 s. Finally, the sections were dehydrated and dried and examined under microscope.

Ethics approval was sought and approved by the Ethics Committee of Fujian Provincial Hospital (Ethics Approval Number K2022-09-103). The patients/participants provided their written in-formed consent to participate in this study.

Development of IGFBP3-Related Risk Score (IGRS)

R language was employed to analyze and acquire the top 200 IGFBP3-related genes (r > 0.2, p < 0.05) in each dataset (Table S1), and then the IGFBP3 and intersection genes were input into the LASSO Cox analysis in the TCGA to construct the prognosis model. The IGRS was calculated follow the formula. ${\displaystyle \text{IGRS}=\sum _i^n\left(\text{Coefficient of}\left(i\right)\times \text{Expression of gene}\left(i\right)\right).}$

Establish and evaluate a nomogram

Based on the R “rms” package, a nomogram incorporating age, gender, tumor grade, tumor stage and IGRS was constructed to predict 1-, 3-, and 5-year survival. Simultaneously, corresponding calibration curves were also plotted to assess the calibration of the nomogram. According to the C-index, the accuracy between nomogram and other prognostic factors was assessed (Chen et al., 2021). Additionally, the decision curve analysis (DCA) was conducted by the “DCA” package to measure the net clinical benefits of various forecasting models (Vickers et al., 2008).

GSEA analysis

DEGs between the low and high IGRS subgroups were identified using “limma” R package, with the standards of |log₂(FC)| > 0.5 and adjusted p < 0.05. The GESA analysis was carried out using the Hallmark and C2 KEGG gene sets v7.4, which were used in conjunction with the GSEA software (version 4.1.0), with p < 0.05 and a FDR of < 0.25 were considered statistically significant.

Immune profile analysis and immunotherapy response analysis

The infiltrating scores of 24 immune cell subtypes was calculated by the IMMUNCELL AI algorithm (Miao et al., 2020). The immune/stromal scores (ImmuneScore and StromalScore) of the LIHC samples were estimated by ESTIMATE algorithm based on given gene expression profile in FPKM or normalized log2 transformed values (Yoshihara et al., 2013). The tumor immune dysfunction and exclusion (TIDE) was calculated to assess the immunotherapy responses in TCGA and validated in the ICGC cohort, as described previously (Jiang et al., 2018).

Statistical analysis

The R version 3.6.3 (R Core Team, 2020) ggplot2 package was used for visualizing the receiver operating characteristics (ROC) curve. The “survival” package was employed to analyze the survival prognosis with the median value of a marker. Independent prognostic factors analysis was conducted by the univariate Cox regression method and multivariate Cox regression method. The Pearson method was used for correlation analysis. p < 0.05 was considered statistically significant.

The expression profiles of IGFBP3 in HCC

To investigate the potential role of IGFBP3 on HCC, we firstly determined the expression profiles of IGFBP3 in HCC sample. The plot indicates that the gene expression of IGFBP3 was relatively downregulated in the HCC samples compared with normal samples (Figs. 2A–2B, all p < 0.001). Moreover, we validated the down-regulation of IGFBP3 in three independent GEO datasets (GSE54236, GSE14520 and GSE76427) (Figs. 2C–2E, all p < 0.001). Furthermore, decreased mRNA expression profile of IGFBP3 was also confirmed in Huh7 and LO2 cell lines (Fig. 2F, p < 0.01). Additionally, HCC tissues showed a significantly decrease of protein expression of IGFBP3 (Fig. 2G, p < 0.001), which was also confirmed (Fig. 2H) by the IHC analysis.

Diagnostic value of IGFBP3 and its relevance to clinical features

A correlation analysis was then performed between IGFBP3 and corresponding clinical characteristics. Statistical significance between age groups was determined using the Wilcoxon rank sum test (p = 0.023). Furthermore, as can be seen in Table 1, the expression level of IGFBP3 showed statistically significant correlation with gender (p = 0.002), T stage (p < 0.001), pathologic stage (p < 0.001), histologic grade (p = 0.005), AFP (p < 0.001) and vascular invasion (p = 0.008). Results also indicated that female patients had higher IGFBP3 levels compared to males (Fig. 3B, p < 0.05), and patients with vascular invasion had higher IGFBP3 expression than those without (Fig. 3H, p < 0.01). In addition, further results revealed significant correlations between the expression level of IGFBP3 and T stages (Fig. 3C), N stages (Fig. 3D), pathologic stages (Fig. 3F), and histologic grades (Fig. 3G) (all, p < 0.05). However, no significant relationships between the expression of IGFBP3 and the age (Fig. 3A), M stage (Fig. 3E) were identified (all, p > 0.05).

To evaluate the diagnostic significance of IGFBP3, ROC analysis was performed based on the TCGA cohort (Discovery cohort). As showed in Fig. 3I, IGFBP3 exhibited powerful diagnostic ability for HCC with an AUC of 0.927 (95% CI [0.902–0.951]). Moreover, the ROC curves of two testing cohorts (GSE14520 and GSE76427) showed that IGFBP3 levels for diagnosing HCC yielded AUCs of 0.934, and 0.910, respectively (Figs. 3J–3K). In addition, in the validation cohort (ICGC-LIRI), IGFBP3 also displayed highly effective in discriminating HCCs from normal samples (Fig. 3L, AUC = 0.912, 95% CI [0.883–0.942]). In order to further analyze its early diagnostic value, we first analyzed the difference of its expression in different stages of the disease. Based on the TCGA cohort, with progressing tumor stages, IGFBP3 gene expression increased (Fig. S1A). According to ROC analysis, IGFBP3 was extremely effective in discriminating early tumor pathologies (stage I and stage II) from normal (Figs. S1B and S1C). These results suggested that IGFBP3 is downregulated in LIHC and can be used as a valuable diagnostic biomarker for LIHC.

Construction and validation of IGFBP3-related risk score

It has been shown that in NSCLC, IGFBP3 mediates growth inhibition and induction of apoptosis to exert a tumor suppressive effect (Hochscheid, Jaques & Wegmann, 2000). Moreover, IGFBP3 has been reported to hinder aggressive growth of HCC in children (Regel et al., 2012). In addition, IGFBP3 has been convinced to correlate with patients response to radiotherapy and chemotherapy in glioblastoma  (Zhao et al., 2011). Given the potential role of IGFBP3, we hypothesized that IGFBP3-related genetic features might be valuable for predicting the prognosis and treatment of hepatocellular carcinoma. Based on Pearson correlation analysis, we first performed analysis in the three datasets (TCGA, GSE14520, and GSE76427), and finally filtered out 10 significantly correlated genes (Fig. 4A). After the LASSO regression analysis, we obtained eight key genes, namely, IGFBP3, RGS2, IER3, PFKFB3, ENO2, FZD1, JUNB, and PELI2 (Figs. 4B–4C). Next, the IGRS was built according to the expression of key genes and their Cox regression coefficients (Table 2). Figure 4D showed the IGRS, survival status, and expression of the eight model genes between high- and low-risk groups in the TCGA dataset. The results of survival analysis suggested that the high IGRS group had a worse survival outcome than the low IGRS group (Fig. 4E; p < 0.001). For overall survival time prediction, IGRS yielded the AUC values of 0.709 at 1 year, 0.705 at 3 years, and 0.716 at 5 years (Fig. 4F). In validating cohorts (GSE14520), the IGRS, survival status, and the expressions of eight model genes were presented in the Fig. 4G, which is similar to the results with TCGA cohort. Survival results also showed a significantly worse survival outcome in the high IGRS group than in the low IGRS group (Fig. 4H; p < 0.001). As can be seen from Fig. 4I, the predicted AUCs of 1, 3, and 5 years were 0.641, 0.682, and 0.592, respectively.

To further define whether IGRS was an independent prognostic factor for OS, we first performed univariate Cox regression analysis. As can be seen from the results, IGRS were significantly associated with OS in both the TCGA [Hazard ratio (HR) = 3.878, 95% CI [2.313–6.502], p < 0.001] and GSE14520 datasets (HR = 3.982, 95% CI = 2.024–7.835, p < 0.001) (Figs. 5A–5B). Afterwards, multivariate Cox regression analysis of both the training and validation cohorts showed that IGRS was an independent prognostic factor for OS (HR = 2.804, 95% CI [1.608–4.890], p < 0.001, and HR = 2.639, 95% CI [1.262–5.519], p = 0.010, respectively).

Establishment of a predictive nomogram

Nomograms are commonly used to quantify risk in individuals. Currently, a nomogram was built according to age, gender, tumor stage, tumor grade, and IGRS (Fig. 6A). The nomogram model with C-index values of 0.684 indicated to have favorable discrimination abilities (Fig. 6B). Additionally, the nomogram showed relatively good agreement with observation in predicting 1-, 3-, and 5-year survival outcomes (Figs. 6C–6E). The DCA curves revealed that the nomogram demonstrated a net benefit over age, sex, grade, stage, and IGRS in terms of 1-, 3-, and 5-year OS (Figs. 6F–6H). In summary, IGRS-based nomogram can predict the short- and long-term OS of HCC patients and help clinical management.

GSEA of IGRS–related signaling pathways

Based on the |log₂FC| ≥ 0.5, FDR < 0.05, the 2021 up-regulated and 437 down-regulated DEGs was identified between the two groups (Fig. 7A). The expression heatmap of the top 60 DEGs was shown in Fig. 7B. The GSEA showed significant enrichment of signatures associated with apoptosis, cell cycle, lysosomes, MAPK signaling pathways, and WNT signaling pathways in the high IGRS group. Moreover, high-risk individuals exhibited enriched expression of the mTOR signaling pathway, JAK STAT signaling pathway, p53 signaling pathway, ERBB signaling pathway, and cancer pathways (Fig. 7C). Interestingly, we found that low-risk was significantly enriched for some metabolism-related pathways, such as, the fatty acid metabolism pathways (Fig. 7D). These results suggested that the two risk groups have different pathway activation states, which may account for the different survival rates.

Immune profile and prediction of treatment style of IGRS-based HCC groups

According to the results of GSEA analysis, immune-related pathways were found to significantly enriched in the high IGRS group (Fig. 8A). In addition, our results showed that the ESTIMATE score and immune score were significantly higher in the high IGRS group compared with the low IGRS group (Figs. 8B–8C).

Recent studies have provided evidence that high expression of checkpoint genes indicated a more sensitive immunotherapy response  (Li, Chan & Chen, 2019). In the current study, we found the elevated expression of CD274, CTLA4, HAVCR2, TIGIT, LAG3, PDCD1in the high IGRS group (Fig. 8D, all p < 0.05). We also observed high IGRS group also expressed higher levels of key chemokines and their receptors (Fig. 8E). A significant increase was also found in MHC-I and MHC-II component levels in the group with high IGRS (Fig. 8F). Detailed differences in immune cell subtypes were further analyzed between the two groups. We found that CD4 naive cells, cytotoxic, exhausted, type 1regulatory T cells (Tr1), natural regulatory T cells (nTregs), iTregs, Th1, Tfh, NK, NK T (NKT), DC, CD4+ T, and CD8+ T cells have a high prevalence in the high-risk group (Fig. 8G, all p < 0.05). In contrast, CD8 naive cells, Th17, monocyte and gamma delta cells were more predominant in the low-risk group (Fig. 8G, all p < 0.05). We then analyzed the correlation between the IGRS and Tumor Immune Dysfunction and Exclusion (TIDE), which are recognized as immunotherapy predictors. We found that patients in the IGRS high group tended to achieve higher TIDE scores (Fig. S1), which proposed that patients with low IGRS scores may benefit from immunotherapy.