Hepatocellular Carcinoma Outcome Is Predicted by Expression of Neuronal Calcium Sensor 1

Hepatocellular carcinoma (HCC) is one of the deadliest cancer types worldwide (1), with rapidly increasing rates in the United States (2, 3). Currently, options to cure or palliate HCC are limited to resection, transplantation, and ablative and transarterial therapies, with therapeutic choices dictated by both tumor burden and the extent of hepatic dysfunction. There is only one FDA-approved first-line, palliative systemic chemotherapeutic agent, an oral tyrosine kinase inhibitor (sorafenib). This treatment is reserved for patients who have failed or progressed after local therapies, offering an overall survival advantage of 3 months with a significant adverse side effect profile (4). To develop and optimize new treatment modalities for HCC, novel biomarkers are needed that can predict tumor aggressiveness and identify druggable targets that can be exploited to create innovative therapies.

Calcium is long known to be a major regulator of cell-cycle progression and cellular proliferation, events promoted by growth factors that often cause oscillatory patterns of Ca²⁺ release from ER stores and Ca²⁺ entry from the extracellular space (5–7). Furthermore, recent clinical and basic research evidence implicates a Ca²⁺-driven pro-metastatic phenotype in aggressive tumors with an unfavorable prognosis (8–10). To date, only one component of Ca²⁺-signaling pathways has been proposed as a potential biomarker for HCC. High expression of S100 Ca²⁺-binding proteins is reported to confer an aggressive phenotype in undifferentiated liver cancers (11–14). Although less well characterized, other members of the Ca²⁺-signaling cascade have been implicated in cancer progression (15).

Neuronal Calcium Sensor 1 (NCS1) is a clinically validated prognostic biomarker in two independent breast cancer cohorts, that is closely related to states of dysregulated Ca²⁺-signaling (10). NCS1, a member of the calmodulin superfamily of EF-hand Ca²⁺-sensing proteins, is involved in various cellular processes, including exocytosis (16), regulation of voltage-gated Ca²⁺-channels, neuroplasticity (17), and calmodulin-mediated signaling (18). NCS1 interacts with inositol 1,4,5-triphosphate receptor type 1 (InsP3R1) and leads to an increase in InsP3-dependent Ca²⁺-release (19), especially in an overexpression scenario (20, 21). NCS1 also is upregulated in injured neurons, acting as a survival factor via activation of PI3K/AKT pathway (22). These functional properties and its expression in different extra-neuronal tissues (23–26) together with its likely role in breast cancer development make NCS1 a possible contributor to carcinogenesis in other tumor types as well.

To investigate the possible role of NCS1 in Ca²⁺-mediated tumor progression and metastasis in HCC, we analyzed The Cancer Genome Atlas (TCGA) Liver Hepatocellular Carcinoma (LIHC) dataset to determine whether there is an association of NCS1 with patient survival and whether NCS1 can distinguish among patients, either by sub-population or etiologic risk factor. In addition, we built a regression model to test for significant coexpression of NCS1 using the approximately 60,000 RNA-sequenced transcripts included in the TCGA LIHC dataset. Using information from the regression model, we identified pathways involved in Ca²⁺ and NCS1 signaling that influence tumor progression and can be pharmacologically manipulated.

The results shown here are in part based upon data generated by the TCGA Research Network: http://cancergenome.nih.gov/. TCGA LIHC, breast cancer (BRCA), and lung cancer (LUAD and LUSC) expression data and clinical annotations were download from the NIH's Genomic Data Commons (GDC) Data Portal (https://portal.gdc.cancer.gov) on October 10, 2017. Plain text files with FPKM values for 60 483 transcripts per patient were joined with corresponding clinical annotations and gene names were added using unique ENSG identifiers (ensembl.org, GRCh38.p11). NCS1 and LIMK1 expression values and clinical annotations related to the International Cancer Genome consortium (ICGC) LIRI-JP liver cancer cohort were downloaded from the ICGC data portal (https://icgcportal.genomics.cn) on February 22, 2018. All data transformation and analyses were performed with R Statistical Programming Language, Version 3.4.0 (2017-04-21; ref. 27). Additional information regarding sample processing, patient characteristics, data acquisition and global mutational and expressional analyses are available with The Cancer Genome Atlas Research Network's publications of these datasets (28–31). Data for NCS1, GAPDH, RPS18, HPRT1 and ACTB expression levels in various tissues that lack HCC (assumed to be healthy for the purposes of this study) were accessed through the websites of The Genotype-Tissue Expression (GTEx) Project (https://www.gtexportal.org/home/) as well as The Human Protein Atlas (HPA; https://www.proteinatlas.org) (both on November 1, 2017). Pre-processed gene expression data of 1,019 cancer cell lines are publicly available through Broad Institute's Cancer Cell Line Encyclopedia (CCLE; https://portals.broadinstitute.org/ccle) and were downloaded on January 10, 2018.

The between-group comparisons of overall survival were performed by means of two-sided log-rank tests stratified according to NCS1 (<75% vs. ≥ 75%), WDR5 (<80% vs. ≥ 80) and LIMK1 (<60% vs. >60%) expression levels. Cutoff values for dichotomizing continuous expression values to “high” and “low” were derived from receiver operating characteristics (ROC) curves using Youden's J statistic. A Cox proportional-hazards model that included dichotomized NCS1, LIMK1, and WDR5 expression as covariates was used to estimate hazard ratios (HRs) and their associated 95% confidence intervals (CIs). The Kaplan–Meier method was used to estimate survival curves and a log-rank test to compare survival between distinct expression subgroups. Survival plots and forest plots were generated using the “survminer” R package (32). In all analyses, patient samples were excluded whenever survival data or expression values were missing (<10% per cohort and analysis).

To assess transcriptional coregulation of the 60,482 transcripts with NCS1 expression in LIHC, BRCA, LUAD, and LUSC datasets, simple linear regression analysis was performed. P values derived from these analyses were corrected using Bonferroni's method and top hits were identified by means of corrected P values. We pre-specified a threshold for Pearson's correlation coefficient of 0.6 or higher to indicate a strong correlation, mainly to reduce the false positive rate (33). Again, a univariate linear regression model was used to compute the relationship between NCS1 and LIMK1 or WDR5 in the TCGA, ICGC, and CCLE datasets.

Exploratory Gene Association Networks (EGAN; ref. 34) software was used to combine NCS1-correlated genes and database-derived gene and pathway annotations into meaningful network plots. A list of 40 genes, comprising the top 10 hits from all 4 datasets, was examined for overlaps with pathways from the Molecular Signature Database (MSigDB, http://software.broadinstitute.org/gsea/msigdb) and by means of EGAN software. Results with an FDR q < 0.05 were called significant and listed as relevant. Welch's two sample t test was used to test for statistical differences of sample means if not otherwise specified and Fisher's exact test was used to analyze contingency tables.

A commercially available HCC microarray with 49 cases (2 cores per case, 41 tumors with matched non-neoplastic adjacent tissue, 6 tumors without matched non-neoplastic tissue, and 2 normal liver tissues without corresponding tumors) was purchased from US Biomax, Inc. (OD-CT-DgLivT10-024, HLiv-HCC180Bch-01). Immunofluorescence staining for NCS1 and arginase was performed and automated quantitative analysis (AQUA) with arginase as a marker for tumor tissue compared with stromal tissue was used to quantify NCS1 protein levels in tumors and adjacent non-neoplastic tissues, as previously described (10). Biochemical indicators (AFP, CA19-9, CEA, hepatitis B/C markers, Creatinine, BUN) and tumor stage data are available on the supplier's website.

To assess physiological NCS1 expression in normal liver tissue from healthy individuals, we explored two databases; GTEx Project and HPA. Analysis of both of these datasets suggests very low overall NCS1 levels in healthy human liver (mean NCS1 expression levels derived from RNA-sequencing specified as 0.9 TPM and 0.7 RPKM in 10 and 119 specimens in the HPA and GTEx datasets, respectively). Furthermore, when compared with other organs, liver expression levels rank last and second to last among 31 and 36 analyzed tissues in the HPA and GTEx datasets, respectively. Because interexperimental comparison of expression units is usually not possible, we calculated the ratio of NCS1 to the housekeeping gene β-actin (ACTB) in HPA, GTEx, and TCGA LIHC datasets (Fig. 1A). The ratio of NCS1 to ACTB was 1.7- and 1.4-fold higher in cancerous tissue (P < 0.01 for the comparison of TCGA vs. GTEx and TCGA vs. HPA, respectively). To validate the choice of ACTB for normalization, three additional housekeeping genes (HPRT1, GAPDH, and RPS1) were used to normalize the NCS1 levels (Supplementary Fig. S1). An increase in NCS1 expression levels in the HCC samples was observed with all housekeeping genes. These results suggest NCS1 levels increase in liver tissue during HCC development.

To further examine this hypothesis, a commercially available microarray composed of HCC samples and adjacent non-neoplastic tissues was analyzed using immunofluorescence to calculate AQUA scores. Matched adjacent non-neoplastic tissues were available for most tumor samples (41/47) on the array. A significant (P < 0.0001) upregulation of NCS1 was observed in tumors (Fig. 1B), with 59% (24/41) of them showing increased NCS1 AQUA scores compared with matched non-neoplastic samples (Fig. 1C). However, in all of the samples, the levels of NCS1 in the non-neoplastic tissue is low and there is a very small range, suggesting that baseline NCS1 will not be a predictor for development of HCC and that high baseline NCS1 expression cannot be a prerequisite for the development of HCC. This is more clearly shown by the relationship between AQUA scores in tumors and non-neoplastic samples (Fig. 1D). Each data point is the expression of NCS1 in tumors and non-neoplastic samples in the same patient and the scatter of the data points show that there is no correlation. This relationship suggests that baseline NCS1 levels do not predict the expression level of NCS1 in the tumor samples and provides further evidence that increased NCS1 expression is associated with tumor development. An analysis of the limited clinical annotations available for the array showed no association of NCS1 AQUA scores in tumor tissue with either tumor stage (Supplementary Fig. S2B), hepatitis B/C status, tumor markers (CEA, CA19-9, AFP) or any of the patient characteristics (age and gender). No additional information on the spatial relationship between tumor samples and non-neoplastic adjacent tissues were available for analysis.

To test whether NCS1 levels correlate with tumor stage and outcome in the TCGA cohort of liver cancer patients, we examined the relationship between NCS1 expression measured in FPKM units and clinical tumor stage for 377 patients. Although advanced tumors did not show increased NCS1 expression compared with early-stage tumors in this cohort (Supplementary Fig. S2A), in the subgroup of Asian patients who died, average NCS1 expression was significantly increased compared with the survivors (t-test corrected P < 0.05 for male Asian patients, Fig. 2A). In addition, Asian patients had a higher mean NCS1 expression than white patients (t test P = 0.056) and a significantly (P < 0.01) smaller proportion of NCS1^low Asian patients died during the course of their disease (21.4% vs. 39.0% for Asian and white patients, respectively, Fig. 2B), irrespective of tumor stage. This comparison suggests a protective effect of low NCS1 expression in NCS1^low Asian patients. Differences in the expression profiles of proteins associated with HCC among races has recently been described (35). With respect to Ca²⁺-signaling pathways, indel mutations in ryanodine receptors (RYRs; intracellular Ca²⁺ channels) were identified in Thai patients in this study (35). Because the number of subjects from races other than Asians and whites constitute only a small part of the overall sample (Supplementary Table S1), they were excluded from further analyses.

In the next analysis, we examined the prognostic value of dichotomized NCS1 expression levels (to NCS1^high and NCS1^low) with regard to patient survival at the last follow-up and over time. To explore the optimum cutoff for dichotomization of continuous FPKM values, we created a receiver operating characteristics (ROC) curve by calculating sensitivity and specificity for a range of different cutoffs and calculated Youden's J statistic to capture the performance of our binary test (Fig. 2C). As a result, we categorized 25% of patients as NCS1^high and 75% as NCS1^low.

When applying this cutoff value to the survival status of Asian patients (i.e., patients being either dead or alive at the time of last follow-up), the estimated odds ratio of 0.26 shows a significant (Fisher's exact test P < 0.001) association of high NCS1 expression levels with an unfavorable disease outcome (Fig. 2E). However, the relationship between NCS1 expression and prognosis is not true for white patients (Fig. 2D), confirming the result of our prior analysis.

We next estimated survival curves via the Kaplan–Meier method and applied a Cox proportional-hazards model with NCS1 as a single covariate to calculate HRs for all patients in the TCGA LIHC cohort (Fig. 3A), as well as race-stratified for the subgroup of Asian (Fig. 3B) and for the subgroup of white patients (Fig. 3C). With this approach, Asians showed a significantly higher probability of death if categorized as NCS1^high (HR, 3.5; P < 0.0001), whereas the time to death in white patients was independent of their NCS1 status (HR, 0.89; P = 0.68). NCS1 status was also significantly associated with worse outcome in the Kaplan–Meier analysis of all patients (HR, 1.7; P = 0.006), although with a hazard ratio closer to 1 than the ratio calculated for Asian patients only.

To validate the strong relationship between high NCS1 expression and survival of Asian patients in the TCGA LIHC cohort, we then examined the International Cancer Genome Consortium (ICGC) LIRI-JP cohort (36) that comprises 231 HCC samples from Japanese patients. After dichotomizing continuous NCS1 expression levels in NCS1^high and NCS^low as described above, we again found that high expression was significantly associated with worse patient survival over time (Fig. 3D; with an unadjusted HR, 1.96; P = 0.039).

Because the exact signaling pathway that is influenced by NCS1 expression in tumor cells remains to be elucidated, we performed simple linear regression analyses of NCS1 with 60 483 RNA-sequenced transcripts from the TCGA LIHC dataset to assess transcriptional coregulation of certain transcripts with NCS1. After correcting for multiple hypothesis testing, a total of 54 genes were significantly (P < 10⁻³⁰ and Pearson's correlation coefficient above the pre-specified threshold of 0.6) and positively correlated with NCS1 expression (Supplementary Fig. S3; Supplementary Table S2). LIM Domain Kinase 1 (LIMK1) was the most significantly coexpressed gene and remarkably, LIMK1 stands out, separate from all the other genes that are associated with NCS1 (Supplementary Fig. S3 and Supplementary Fig. S4A, NCS1 vs. LIMK1 Pearson's r = 0.76, P < 0.0001). We further validated the strong coexpression by exploring the relationship of NCS1 and LIMK1 in the ICGC LIRI-JP cohort (Supplementary Fig. S4B, Pearson's r = 0.72, P < 0.0001).

With this information, we constructed a coexpression network (37) using EGAN (Fig. 4A; ref. 34) software. We found cytoskeleton organization to be the dominant pathway in the derived network. To further minimize bias and improve robustness of the approach, we analyzed three additional TCGA cohorts (LUAD, LUSC, and BRCA; Supplementary Table S2 for top hits and Supplementary Fig. S5 for gene association networks) using the same methods and included only the 10 most significant results from each cohort in our pathway analysis.

In this expanded cohort, two genes occurred among the top 10 NCS1-associated genes in at least two of four samples, namely LIMK1 (overlapped between two datasets) and WD Repeat Domain 5 (WDR5, overlapped between three datasets). We calculated a Cox proportional-hazards model that includes dichotomized NCS1, LIMK1, and WDR5 expression as covariates and adjusts for patient age and gender to investigate the prognostic significance of these three NCS1 associated genes in the TCGA LIHC cohort. Although not all of the calculated hazard ratios become significant due to the relatively small number of events, all analyses predicted an increased risk of death in patients with high expression of the respective genes compared with patients with low expression (Fig. 4B; Supplementary Fig. S6A shows the respective forest plot for white patients). These effects were stable and independent of patient age or gender. In addition, we examined the genes that overlap with pathways from the Molecular Signature Database (MSigDB). Again, we found that the genes identified in our analysis are related to either Ca²⁺-dependent signaling pathways (cell-cycle progression, cytoskeleton organization), neurogenesis or cancer-related pathways (Supplementary Fig. S6B).

We also examined Broad Institute's Cancer Cell Line Encyclopedia (CCLE) dataset that comprises RNA sequencing derived expression data of 1,019 cell lines to further test and potentially validate the positive correlation between NCS1 and LIMK1 or WDR5. Consistent with our pathway analysis, the relationship between NCS1 and LIMK1 is observed in breast and liver cancers, but not in cell lines derived from hematological malignancies that do not require increased motility and cytoskeleton remodeling for metastasis. In contrast, NCS1 and WDR5 did not show a significant correlation in the CCLE dataset, independent of cancer type (Supplementary Fig. S7A and S7B).

In this study, we systematically examined the prognostic role of NCS1 and transcriptionally associated genes in a liver cancer cohort from a publicly available database. Interestingly, healthy liver tissue exhibited lower overall NCS1 levels compared with HCC, suggesting an upregulation of NCS1 during tumorigenesis. We were able to validate this finding by using a tissue microarray where 59% of tumors harbored NCS1 levels (as assessed by a quantitative immunofluorescence technique) that were increased compared with the matched non-neoplastic tissues. We did not see a correlation between NCS1 levels in these matched adjacent non-neoplastic tissues and NCS1 levels in the respective tumor samples. This result supports our hypothesis that NCS1 is upregulated during carcinogenesis, in a manner independent of basal levels in the healthy organ or tumor-adjacent parenchyma. We suggest that NCS1 expression is a fixed, tumor dependent phenomenon that is associated with risk factors or concomitant liver diseases, irrespective of stage at diagnosis.

By correlating RNA-sequencing–based expression to clinical outcome measures, specifically survival status of patients at the end of follow-up and time to death, we demonstrated a significantly worse prognosis for Asian patients with high NCS1 expression compared with Asians with low NCS1 expression. We validated this finding by investigating an additional, independent cohort of Asian HCC patients. Although we found NCS1 to be predictive of survival in an analysis of all TCGA samples, we hypothesize that the observed effect is mainly driven by the subpopulation of Asian patients because we did not see an association of NCS1 with survival for white patients. In addition to these analyses, we identified an associated signaling network whose central components further improved survival predictions. The expression pattern of these genes suggests a path for the development of prognostic biomarkers and potential targeted therapies.

One possible cause for the observed race-specific differences regarding survival might be the overall higher expression of NCS1 among Asian cancer patients. This differential expression may be caused by hepatocellular carcinoma-associated risk factors, including hepatitis or alcohol exposure, although a mechanistic explanation is beyond the scope of this study. Nevertheless, it is known that HCCs differ between Asian and white patients in terms of genomic properties (35) and with respect to pro-carcinogenic risk factors (38). In this regard, our study adds to an increasing body of literature investigating possible contributors to these phenomena. Importantly, our study suggests a specific cellular mechanism, dysregulated Ca²⁺ signaling, as a hallmark of HCC in a subgroup of patients.

To facilitate a better understanding of the biology of tumors expressing high levels of NCS1, we investigated additional cancer datasets, two lung cancer and one breast cancer datasets, for significant coregulation of NCS1 with other transcripts. We found a range of pro-oncogenic genes to be highly expressed in NCS1^high patients and identified the associated pathways to be cytoskeleton organization, neurogenesis and cell-cycle progression. Furthermore, most of the identified genes are involved in cancer- and metastasis-related pathways, underscoring the validity of the analysis.

Two genes appeared to be NCS1-associated in at least two of the four examined TCGA datasets (WDR5 and LIMK1). WDR5 is a protein shown to be Ca²⁺ sensitive in amphibians (39) and it has a critical role in recruiting the highly oncogenic transcription factor MYC (40) and epigenetic modifiers MLL1–4 (41, 42) to chromatin (43). This study is the first time WDR5 has been proposed to be associated with NCS1 as a crucial (dys)regulator of intracellular Ca²⁺-signaling. The WDR5/MLL complex has been suggested as a contributor to an unfavorable disease outcome in HCC (44) and MLL4 is a known target of HBV integration, especially in Chinese patients (45). However, WDR5 ranks only among the top 500 genes when analyzing the TCGA LIHC dataset. To include WDR5 as a critical component of tumorigenesis would require a robust validation of a functional role of WDR5 beyond mere association, especially as it is located in close chromosomal proximity to NCS1 and was not correlated in an analysis of RNA sequencing data from the CCLE.

LIMK1 is a promising therapeutic target that was previously studied with regard to metastasis and appeared prominently in our analyses of TCGA and CCLE datasets. Cellular motility, a hallmark of metastatic cancers, is mainly facilitated via actin polymerization (46) and LIMK1 is a key regulator of this process (47). In addition, Ca²⁺-dependent LIMK1 activation was observed in neuronal outgrowth (48). The second most significantly coexpressed gene in the analyzed liver cancer cohort is ARHGEF10, a RhoA guanine nucleotide exchange factor that acts upstream of LIMK1 (49). From the potential functional combination of these proteins along with components of our coexpression network (Fig. 4A), we hypothesize an increased activation of the RhoA—Rho-associated protein kinase (ROCK)—LIMK1 signaling axis in tumor cells exhibiting increased NCS1 levels (50). This functional complex confers motility and allows tumor cells to interact with the stroma and microenvironment, and expressional upregulation would ultimately lead to decreased patient survival. Because potent inhibitors of LIMK1 have been developed that can reduce cellular motility in vitro (51), LIMK1 is a potential clinical target in advanced HCC, possibly in combination with inhibitors that address additional components of a larger signaling network. Furthermore, our results suggest that the NCS1/LIMK1 signaling network can be modeled using cancer cell lines.

Other Ca²⁺-dependent pathways that have been causally related to HCC development (e.g., the lysophosphatidic acid pathway; ref. 52) may be directly or indirectly associated with NCS1 expression. However, these pathways were not identified in the cohort included in this study.

The main limitations of this study are the relatively small number of Asian patients in the TCGA LIHC cohort, the focus on pre-processed RNA-sequencing as the primary measure of protein expression, and the limited availability of patients' past medical histories. Our hypothesis that NCS1 increases during hepatocarcinogenesis is based on a TMA with HCC samples and non-neoplastic adjacent tissues. To test this hypothesis, longitudinal studies are required to explore the role of NCS1 during the evolution of HCC in greater depth. This is especially import with regard to recent evidence suggesting that liver tissue adjacent to a primary liver tumor might have expression patterns or phenotypes different from normal hepatocytes (53).

Moreover, the excellent overall survival rate of about 65% after 10 years in the TCGA study suggests a bias in this cohort toward patients that were treated with a curative intention and thus had on average a longer survival than usually observed with HCC. Nevertheless, our study points to cytoskeleton remodeling, possibly facilitated by an interaction of NCS1 and the RhoA/ROCK/LIMK1 pathway, as an important determinant of aggressive HCC. The question whether disruption of this pathway applies only to a distinct subgroup of patients or to all HCC patients with NCS1 levels beyond a certain cutoff value, will be addressed in subsequent studies.

In summary, using a database-focused proof-of-concept study to analyze the predictive capacity of the multifunctional Ca²⁺-binding protein NCS1 in liver cancer, we found that high expression of NCS1 and transcriptionally coregulated proteins is associated with an unfavorable clinical outcome in HCC. This study lays the foundation for utilizing NCS1 as a prognostic biomarker in prospective cohorts of HCC patients. Furthermore, our work can serve as a starting point to assess the identified signaling complex and demonstrate a functional contribution to tumorigenesis and metastasis.

B.E. Ehrlich has ownership interest in and is a consultant/advisory board member for Osmol Therapeutics. No potential conflicts of interest were disclosed by the other authors.

Conception and design: D. Schuette, B.E. Ehrlich

Development of methodology: D. Schuette

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): L.M. Moore

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Schuette, L.M. Moore, B.E. Ehrlich

Writing, review, and/or revision of the manuscript: D. Schuette, M.E. Robert, T.H. Taddei, B.E. Ehrlich

Study supervision: B.E. Ehrlich

D. Schuette received a scholarship from the German National Merit Foundation. NIH grant 5P01DK057751 (B.E. Ehrlich and M.E. Robert) supported this work.

Helpful discussions with Michael H. Nathanson (Yale University) and Joseph K. Lim (Yale University) are acknowledged.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.