Dynamic network biomarker indicates pulmonary metastasis at the tipping point of hepatocellular carcinoma

Developing predictive biomarkers that can detect the tipping point before metastasis of hepatocellular carcinoma (HCC), is critical to prevent further irreversible deterioration. To discover such early-warning signals or biomarkers of pulmonary metastasis in HCC, we analyse time-series gene expression data in spontaneous pulmonary metastasis mice HCCLM3-RFP model with our dynamic network biomarker (DNB) method, and identify CALML3 as a core DNB member. All experimental results of gain-of-function and loss-of-function studies show that CALML3 could indicate metastasis initiation and act as a suppressor of metastasis. We also reveal the biological role of CALML3 in metastasis initiation at a network level, including proximal regulation and cascading influences in dysfunctional pathways. Our further experiments and clinical samples show that DNB with CALML3 reduced pulmonary metastasis in liver cancer. Actually, loss of CALML3 predicts shorter overall and relapse-free survival in postoperative HCC patients, thus providing a prognostic biomarker and therapy target in HCC.


Reviewer's Comment 1:
It is really annoying that the actual transcriptomic data (so called "whole genome profile") on which the whole study is based is not explained at all. Is this microarray? RNASeq? How are raw data processed?
This has to be clearly explained.

Authors' Response:
Thanks for your comment. To address this question, we added the detailed description of transcriptomic data in Methods of the revised manuscript, as follows.
"Whole-genome expression profile To assess whole-genome expression, we individually extracted total RNAs from liver tumours of orthotopic xenograft mice using TRIZOL Reagent (Life technologies) and checked for a RIN number to inspect RNA integrity by an Agilent Bioanalyzer 2100 (Agilent technologies). The qualified total RNAs were further purified by RNeasy micro kit (QIAGEN) and RNase-Free DNase Set (QIAGEN). Then, total RNAs were amplified, labelled and purified by using GeneChip 3'IVT Express Kit to obtain biotin labelled cRNA (Affymetrix). Array hybridization and wash were performed with constant rotation on the PrimeView Human Gene Expression Assay (Affymetrix). Microarrays were scanned by GeneChip Scanner 3000 (Affymetrix) and Command Console Software 3.1 (Affymetrix) with default settings. Raw data were normalized by robust multiarray analysis (RMA) algorithm, Gene Spring Software 11.0 (Agilent technologies). All microarray data were deposited in the Gene Expression Omnibus database (http://www.ncbi.nlm.nih.gov/geo/) Accession Number GSE94016.
To identify differentially expressed genes, we compared gene expression intensities among samples at two different time points using Welch's t-test with two-tailed p value < 0.05 which was adjusted by False Discovery Rate (FDR) for multiple testing."

Reviewer's Comment 2:
Overall, the manuscript is kind of sloppy. For example, "Pearson's correlation correlation" in Supplementary Fig. 1 or "to overexpress CALML3 in HepG2 cells, with low metastatic potential, and to decrease expression in HCCLM3 cells, with high metastatic potential "(lanes 287-288): this is the reverse! Line 368: "CALML3 expression and so on...".

Authors' Response:
We appreciate that you checked our paper carefully and helped us to improve the quality of this paper.
According to your comment, we double-checked and carefully revised our paper. In particular, we made the following revisions in the revised manuscripts: 1. "Pearson's correlation correlation" was replaced with "Pearson correlation coefficient" in Supplementary Fig. 1. Meanwhile, to address your next question, we added the description of DNB algorithm given in Supplementary Fig. 1, which is explained in the response of the next comment.
2. We made inappropriate descriptions in the previous version of the manuscript, which have been corrected in the revised manuscript as follows.
"In addition, to further verify the role of CALML3 as a tumour suppressor, we overexpressed CALML3 in HCCLM3 cells with high metastatic potential and decreased CALML3 expression in HepG2 cells with low metastatic potential by lentiviral-mediated expression and RNA interference, respectively (Supplementary 3. We re-described the result for easy reading, i.e., "In total, sixteen Clinical characteristics (e.g. tumour encapsulation, vascular invasion, tumour number, tumour size, BCLC stage, intratumoural CALML3 expression and so on) were analysed to identify factors associated with overall survival (OS) and relapse-free survival (RFS)."

Reviewer's Comment 3:
The "detailed" algorithm of the DNB method is not detailed enough. What tests are used for "significantly high deviations"? What tests are used to test for deciding whether a correlation is "high compared to control"? What would be the control here?

Authors' Response:
According to your suggestion, we added the details of DNB algorithm in Supplementary Fig. 1. 1. What tests are used for "significantly high deviations"?
For deviation test, the criterion is for one gene more than 2 times of the ratio (SDd) of standard deviation at t to at t0.
2. What tests are used to test for deciding whether a correlation is "high compared to control"? What would be the control here?
We considered gene expression data of samples at the second week after orthotopic implantation (i.e. t=t0=0) as control (see Supplementary Fig. 1). Here, "compared to control" means that all three criteria

Reviewer's Comment 4:
Lane 218: the authors write: "This result is consistent with previous results analysing the dynamics of gene expression" What are the (unreferenced) "previous results" here?

Authors' Response:
Here, we revised this following sentence with clear illustration in light of the comment.
"This result was consistent with the morphological alterations in orthotopic xenograft mice (Fig. 1b, c) and our assumption from the dynamics of gene expression (Fig. 2b)."

Reviewer's Comment 5:
How many genes are finally being considered as part of the DNB? A list of those genes could be useful (if not too long)?

Authors' Response:
The DNB including 334 genes was detected by our method and it is too many to be listed in the main text.
Thus, we added the following description and also Supplementary Table 2 in the revised manuscript.
"Meanwhile, we obtained the corresponding leading group as DNB which was composed of 334 genes or proteins (Supplementary Table 2)."

Reviewer's Comment 6:
The procedure for choosing CALML3 among the DNB genes is really not properly described. Are not DNB genes a subset of DEGs? If so (and this seems to be the case from Figure 3a) why is it being taken into account? Figure 3a is confusing. What percent is represented among the y axis? This should be made clear in the legend. Why is the PIK3R5 gene not considered as a better candidate that the CALML3 gene?

Authors' Response:
1. The procedure for choosing CALML3 among the DNB genes is really not properly described.
There are 334 genes in the DNB group. It is really a difficult task to measure the importance of a gene in a molecular network due to so many involved factors, and there may be many different ways to do it. To validate the important role of DNB, in this manuscript we considered the following criteria with the priority to rank these DNB genes for conducting biological experiment and added the detailed description in Methods.

"Ranking scheme for DNB members
We ranked the DNB genes according to the following criteria with four priorities.
(1) Priority one (the first bar in Fig. 3a). Rank DNB genes for their importance at the network level. Hub genes are considered to play a central role in the molecular network (including DNB genes and DEGs) during metastasis initiation. A hub from DNB genes is a gene with a number of links (i.e. degree) that considerably exceeds the average in the network. Thus, we mapped both DNB genes and DEGs into the molecular network and then individually counted the total number of DEGs directly linked with each DNB gene, i.e., obtained the degree of each DNB gene. This criterion is represented as a percentage or ratio of DEGs to its neighbouring genes.
(2) Priority two (the second bar in Fig.3a). Rank DNB genes based on their functional roles at biological pathway level. We counted the number of metastasis-associated pathways involving each DNB gene.
Thus, we ranked these DNB genes according to the total number of their involved KEGG-annotated pathways. In this work, 163 pathways were considered to be related to metastasis. These 163 pathways include most metabolism processes (e.g., carbohydrate, energy, lipid, amino acid), and all genetic information processing, cellular processes, and signalling pathways (Supplementary Table 3).
(3) Priority three (the third bar in Fig.3a). This criterion indicates if or not a DNB gene belongs to one of DEGs, i.e., 1 (yes) or 0 (no).
(4) Priority four (the fourth bar in Fig.3a). This criterion indicates if or not a DNB gene belongs to one of the six clusters, i.e., 1 (yes) or 0 (no). If yes, this criterion further indicates which cluster this DNB gene belongs to. In other words, this criterion describes how each DNB gene is related to the role of suppressors or oncogenes at metastasis initiation. Specifically, we checked which DNB genes overlapped with DEGs in Clusters 1 to 6 (Fig. 2b). According to the change in patterns at the gene expression level during W2-W4, we considered that the DNB genes in Clusters 1, 3 and 5 might act as (or be related to) suppressors, preventing from metastasis, because they were downregulated before the metastasis; whereas other DNB genes in Clusters 2, 4 and 6 might function as (or be related to) oncogenes, promoting metastasis, because they were upregulated before metastasis." In Results, we also made the corresponding revision, as follows.
"To further understand the underlying roles of DNB genes in metastasis initiation for primary tumour cells, we ranked comprehensively the DNB genes according to the criteria with four priorities, i.e. importance in the networks, importance in functional pathways, differential patterns, and dynamic patterns (see Methods), and selected CALML3 as the top one for further functional study (Fig. 3a, Supplementary Table   2). Certainly, based on the different considerations, there may be other ways to rank the genes." 2. Are not DNB genes a subset of DEGs? If so (and this seems to be the case from Figure 3a) why is it being taken into account? DNB genes were identified only based on the three criteria (i.e., drastic fluctuation at genes expression, increasing correlations between genes within DNB group, and decreasing correlations between genes within and without DNB group). Thus, there is no requirement on differential expressions for those DNB genes. In other words, a DNB gene may be a DEG or may be not a DEG gene. Hence, DNB genes are not necessary to be a subset of DEGs.
However, as mentioned above, we could estimate roughly what roles DNB genes play according to their changes in patterns and how DNB genes are associated with the DEGs at the gene expression level.
Hence, we took this criterion into account. In light of the comment, we added the explanation in the revised manuscript.
3. What percent is represented among the y axis? This should be made clear in the legend.
As shown in Fig. 3a, 4 different bars represent criteria or attributes mentioned above to rank DNB genes.
For the first bar, we used a percentage to represent how many DEGs were associated with the corresponding DNB gene (degree of each gene in the network) against the number of all neighboring genes for this DNB gene. Note that, in order to avoid bias, we only considered all neighboring genes (including non-DEGs) detected in this whole-genome expression profile. Thus, we re-calculated the corresponding data in Supplementary Table S2 as well as re-ranked the DNB genes and updated Fig. 3a.
The second bar shows the total number of DNB-involved pathways. The third and the fourth bars represent individually if or not a DNB gene belongs to one of DEGs and one of the six clusters which changed during 2 to 4 weeks after orthotopic implantation. In addition, we also labelled the corresponding cluster of DEGs to clearly show how DNB gene expression changes during the critical period of metastasis. As mentioned above, there are different scales in this figure. Thus, to avoid the confusion, we removed the y axis and directly labelled information on the corresponding bars as follows (Fig. 3a). The further explanations were added in the revised manuscript. 4. Why is the PIK3R5 gene not considered as a better candidate that the CALML3 gene?
Indeed, DNB as a group, its members may play important roles during metastasis initiation. Here, we intended to choose not only important and but also novel one for further functional study based on our heuristic criteria. Undoubtedly, FGFR4, PIK3R5 and PDGFRA have been reported notably associated with carcinogenesis, which in turn, further validated the significance and effectiveness of our method to identify the important genes. However, there is less report about relations between CALML3 and cancer, especially in pulmonary metastases of hepatocellular carcinoma. Thus, we chose CALML3 to uncover its dysfunctions during pulmonary metastases of hepatocellular carcinoma. Other genes will be studied as our future work. In light of your comment, we added the related explanations in the revised manuscript.

Reviewer's Comment 7:
A more deep question is: the author assume that CALML3 is kind of the driver ("regulating these reversing DEGs at a molecular network level."). Two questions are immediately obvious: how (i.e. by what molecular mechanism?) and how direct? (i.e. not comparing two established cell lines but a time dependent increase/decrease using an inducible system).

Authors' Response:
As suggested by the reviewer, CALML3 may not be a kind of driver although DNB genes are considered to be strongly related to drivers. But rather than drivers, the main purpose of our study is to identify the tipping point of HCC metastasis by dynamic network biomarker that is able to signal the emergence of the disease transition. Thus, based on your comment as well as the suggestion from the reviewer #2, we changed the description of CALML3 from "driver" to "indicator" of pulmonary metastases in the revised manuscript because it indicates the imminent transition of pulmonary metastases based on our data. In other words, we reconsidered CALML3 as an indicator of pulmonary metastases, which is also a suppressor based on the result of present study.
In addition, to further explore the role of CALML3 in HCC metastasis, we performed gain-of-function and loss-of-function studies to detect the effect of CALML3 on oncogenic behaviours of hepatoma cells including cell growth, migration and invasion. Considering that MHCC97L and HCCLM3 (metastatic potential: MHCC97L < HCCLM3) were established from the same parent HCC cell line and showed relative high and low CALML3 expression respectively, we established a doxycycline-inducible CALML3 overexpression system in the HCC cell line HCCLM3 and a CRISPR/Cas9 mediated knock out in the HCC cell line MHCC97L. Treatment with 50 ng/ml doxycycline induced a time-dependent expression of CALML3 protein in HCCLM3/CALML3 cells, detectable as early as 12 h (Fig. 3c). Two CALML3-targeting gRNAs sequences were designed to avoid non-specific effects of the CRISPR/Cas9 system.
Cancer is characterized by sustaining proliferative signalling, evading growth suppressors, resisting cell death, enabling replicative immortality, inducing angiogenesis, and activating invasion and metastasis 38 .
First we examined cell proliferation ability of CALML3 overexpression or knockout cells and the corresponding controls. The result of cell proliferation assay showed that inducible CALML3 expression by Tet-on system significantly inhibited cell growth in HCCLM3 cells, while knockout of endogenous CALML3 by CRISPR/Cas9 technology could promote cell proliferation in MHCC97L cells (Fig. 4a). The results of cell migration and invasion assay showed that compared to the controls, CALML3 overexpression cells displayed markedly decreased ability of migration and invasion, however, CALML3 knockout had an opposite effect on metastasis of HCC cells (Fig. 4b-e). In accordance with assay in vitro, the result of xenografts models established by inducible CALML3 overexpression system in HCCLM3 cells verified that overexpression of CALML3 significantly diminished tumorigenic capacity and pulmonary metastasis. Inducible CALML3 expression group exhibited less weight and volume of tumours, as well as less pulmonary metastasis nodules and incidence than the control (Fig. 4f-i). In addition, to further verify the role of CALML3 as a tumour suppressor, we overexpressed CALML3 in HCCLM3 cells with high metastatic potential and decreased CALML3 expression in HepG2 cells with low metastatic potential by lentiviral-mediated expression and RNA interference, respectively ( Supplementary Fig. 4a). The results also showed that overexpression of CALML3 remarkably inhibited cell proliferation ( Supplementary Fig.   4b), cell migration ( Supplementary Fig. 4c,d), and cell invasion ( Supplementary Fig. 4e,f), whereas downregulation of CALML3 had the opposite effect. We also created xenografts models of the control HCCLM3 and CALML3 overexpression HCCLM3 cells with fluorescent proteins as above (Fig. 1a); CALML3 overexpression significantly inhibited tumour growth and tumour pulmonary metastasis ( Supplementary Fig. 4g,h). Together, the results of these functional assays in vitro and in vivo, suggested mice implanted with HCCLM3 cells (the control and CALML3 overexpression cells). CALML3 overexpression showed less tumour weight. (h) Lung tissues from mice implanted with HCCLM3 cells (the control and CALML3 overexpression) orthotopic transplantation model were examined by fluorescence microscopy (Left panel). CALML3 overexpression HCCLM3 cells showed less lung metastatic nodules (Middle panel) and reduced incidence of lung metastasis (Right panel). *p < 0.05, **p < 0.01.

Reviewer's Comment 8:
In the discussion, the authors write that their approach is "model free". This I think is a misleading assumption: they assume the existence of a dynamical system in which genes-to-genes interaction do exists and can be captured using a simple correlation coefficient. This is a strong assumption which characterizes a model in my view (although I fully agree that they do not try to infer the parameter values of such a dynamical system, they assume its existence).

Authors' Response:
Here, the DNB method is considered as "model free" because the three conditions of DNB were derived from the generic property of nonlinear dynamical systems. It can be proven that data generated from a nonlinear system (defined by ordinary differential equations) should satisfy the three conditions of DNB theory near the tipping point when the system state approaches the tipping point starting from a stable equilibrium. Thus, it is not necessary to infer the unknown parameters. However, as indicated by the reviewer, there may be other complicated issues due to the complexity of biological systems, and thus we removed the "model free" in the revised manuscript.
Proofing a tumour suppressor requires endogenous genetic tests on "loss of function" and "gain of function". Available data could only support a putative tumour suppressor at best. Similarly for metastatic suppression.

Authors' Response:
We appreciate the reviewer's valuable comment. According to your suggestion, we established a doxycycline-inducible CALML3 overexpression system in the HCC cell line HCCLM3 for gain-of-function assays and a CRISPR/Cas9 mediated knock out in the HCC cell line MHCC97L for loss-of-function assays. Those experiments both in vitro and in vivo clearly validated our conclusions of this work (Figs. Supplementary Figs. 3-4). In addition, our study also verified that CALML3 was a suppressor gene in HCC tumourigenesis and metastasis. In light of your suggestion, we have significantly improved our work in the revised manuscript, and the detail revisions in the manuscript are as follows.

3-4 and
"To further explore the role of CALML3 in HCC metastasis, we performed gain-of-function and loss-of-function studies to detect the effect of CALML3 on oncogenic behaviours of hepatoma cells including cell growth, migration and invasion. Considering that MHCC97L and HCCLM3 (metastatic potential: MHCC97L < HCCLM3) were established from the same parent HCC cell line and showed relative high and low CALML3 expression respectively, we established a doxycycline-inducible CALML3 overexpression system in the HCC cell line HCCLM3 and a CRISPR/Cas9 mediated knock out in the HCC cell line MHCC97L. Treatment with 50 ng/ml doxycycline induced a time-dependent expression of CALML3 protein in HCCLM3/CALML3 cells, detectable as early as 12 h (Fig. 3c). Two CALML3-targeting gRNAs sequences were designed to avoid non-specific effects of the CRISPR/Cas9 system. Confirmation of the genotype of two CALML3 -/cell lines (CALML3 -/-#1 and CALML3 -/-#2) was testified by genomic DNA sequencing and Western blot ( Supplementary Fig. 3 and Fig. 3d).
Cancer is characterized by sustaining proliferative signalling, evading growth suppressors, resisting cell death, enabling replicative immortality, inducing angiogenesis, and activating invasion and metastasis 38 .
First we examined cell proliferation ability of CALML3 overexpression or knockout cells and the corresponding controls. The result of cell proliferation assay showed that inducible CALML3 expression by Tet-on system significantly inhibited cell growth in HCCLM3 cells, while knockout of endogenous CALML3 by CRISPR/Cas9 technology could promote cell proliferation in MHCC97L cells (Fig. 4a). The results of cell migration and invasion assay showed that compared to the controls, CALML3 overexpression cells displayed markedly decreased ability of migration and invasion, however, CALML3 knockout had an opposite effect on metastasis of HCC cells (Fig. 4b-e). In accordance with assay in vitro, the result of xenografts models established by inducible CALML3 overexpression system in HCCLM3 cells verified that overexpression of CALML3 significantly diminished tumorigenic capacity and pulmonary metastasis. Inducible CALML3 expression group exhibited less weight and volume of tumours, as well as less pulmonary metastasis nodules and incidence than the control (Fig. 4f-i). In addition, to further verify the role of CALML3 as a tumour suppressor, we overexpressed CALML3 in HCCLM3 cells with high metastatic potential and decreased CALML3 expression in HepG2 cells with low metastatic potential by lentiviral-mediated expression and RNA interference, respectively ( Supplementary Fig. 4a). The results also showed that overexpression of CALML3 remarkably inhibited cell proliferation ( Supplementary Fig.   4b), cell migration ( Supplementary Fig. 4c,d), and cell invasion ( Supplementary Fig. 4e,f), whereas downregulation of CALML3 had the opposite effect. We also created xenografts models of the control HCCLM3 and CALML3 overexpression HCCLM3 cells with fluorescent proteins as above (Fig. 1a); CALML3 overexpression significantly inhibited tumour growth and tumour pulmonary metastasis ( Supplementary Fig. 4g,h). Together, the results of these functional assays in vitro and in vivo, suggested that CALML3 was a suppressor gene in HCC tumorigenesis and metastasis."