Genomic features reveal the clonal evolution in human hepatocellular carcinoma cell lines MHCC97 with stepwise  metastatic  potentials

doi:10.21203/rs.3.rs-32352/v1

Download PDF

Research article

Genomic features reveal the clonal evolution in human hepatocellular carcinoma cell lines MHCC97 with stepwise metastatic potentials

https://doi.org/10.21203/rs.3.rs-32352/v1

This work is licensed under a CC BY 4.0 License

Version 1

posted

You are reading this latest preprint version

Beckground：The series of human Metastatic Hepatocellular Carcinoma (MHCC) cell lines MHCC97 with a continuously enlarged metastatic capabilities have been widely used as a metastatic model in cancer research. However, their genomic features have not been thoroughly investigated so far. This study aims to reveal the genetic events that could drive metastatic clonal evolution through uncovering their genomic aberrations, together with methylomic and transcriptomic alteration.

Results： We investigated the karyotypes on MHCC97L and HCCLM6 HCC cell lines with lower and the highest metastatic capability respectively, and interrogated their complete repertoire of genomic features by performing whole genome sequencing (WGS), whole exome sequencing (WES) and single-nucleotide polymorphism (SNP) genotyping. Mutations including single nucleotide variant (SNV), copy number variation (CNV) and structural variant (SV) were called. In addition, the methylome profiles were analyzed via Illumina Infinium HumanMethylation450 Bead Chip and the transcriptome profiles were analyzed via RNA sequencing (RNA-Seq). The karyotypic examination revealed that MHCC97L cells are more heterogenous than the HCCLM6 cells. The clonal analysis on the two cell lines revealed that a minor subclone with the mutations predisposed to metastatic potential within the MHCC97L become the dominant one within the HCCLM6 cell population. The amplification of CAV1, CAV2 and MET, together with STAT1 deep deletion and DELC1 promoter hypermethylation could be the earlier key drivers occurred in the MHCC97L cells. Subsequently, extra CCL4, CCL26 and TWIST2 amplification, along with the c.1061T>G in SMAD5 and c.1373A>G in RNF169 mutations and deregulatory promoter methylation in the HCCLM6 cells, aggravate the metastatic capability.

Conclusions： Our results revealed that HCCLM6 cells might be derived from a minor subclone within the MHCC97L cell population, where the genetic events, including MET, CAV1 and CAV2 amplification, together with the extra mutations, amplification and deregulatory promoter methylation of some genes, drive the clonal expansion and evolution of those cells with the metastatic capability.

Cancer Biology

Oncology

MHCC97L

HCCLM6

SNV

CNV

Hepatocellular carcinoma (HCC) is one of the most common human malignancies [1]. In these countries with epidemic hepatitis B virus (HBV) infection, HCC is prevalent and leads to a high mortality rate every year. In addition, other risk factors, including aflatoxin contamination, alcoholism, diabetes and hepatitis C virus (HCV) infection, have made HCC being a worldwide health issue [2-4].

As known, metastasis is the leading cause of cancer-related mortality. During the metastatic process, HCC cells, like most other types of cancers, are continuously evolving in a stepwise or punctuated style under selection pressure [5-9]. They may acquire novel genotypes and subsequent phenotypes to migrate to distant territories. Although people have investigated a lot of HCC patients with metastasis, the somatic genetic aberrations driving the cancer metastasis and malignant clonal evolution remain obscure due to the complexity of in vivo tumor environment. Thus, a suitable in vitro metastasis model system is needed to investigate the genetic mechanisms driving HCC metastasis.

Fortunately, the Liver Cancer Institute of Fudan University in China established a continuum of Metastatic Hepatocellular Carcinoma (MHCC)-related HCC cell lines, including MHCC97L, MHCC97H and HCCLM6, derived from an intrahepatic metastatic lesion of a male Chinese HCC patient with HBV infection background [9-13]. In these HCC cell lines, both MHCC97L and MHCC97H cell lines with different metastatic potentials were derived from their parental metastatic cell line MHCC97, where the MHCC97L was separated as a relatively lowly metastatic fraction, compared to its counterpart MHCC97H. The HCCLM6 was produced in the procedures of repeating six rounds of MHCC97H orthotopic implantation, lung metastasis dissection, and subcutaneous implantation. With iterative selection, the HCCLM6 became metastatic even under subcutaneous implantation or footpad implantation [9-13]. Currently, these MHCC97 derivers, MHCC97L, MHCC97H and HCCLM6 with a continuously enlarged metastatic capabilities, have been widely adopted in HCC studies, especial HCC metastasis. Although these cell lines were derived from the same tumor, the heterogeneous cell population within primary HCC tumor may exhibit distinct metastatic potentials. It is possible that these MHCC97 derivers through in vivo gradually selection strategy could partially represent the evolutionary process and different subclones of diversiform HCC cell populations, although the in vivo selection is artificial.

However, it is intriguing how these cell lines with the analogy in their genomic makeups diversified their viabilities in the course of intruding distant territories. We proposed that certain genetic events responsible for the metastatic clonal evolution are involved in these MHCC97 derivers. The elucidation of these genetic events that could drive the metastatic clonal evolution will be helpful to uncover the mechanisms related to HCC metastasis. To explore the genetic events involved in the MHCC97 derivers, here we figure out their genomic aberrations, together with methylomic and transcriptomic alteration, of both MHCC97L and HCCLM6, representing the two terminals of the MHCC97 cell lines. There data will lay a foundation for the understanding of metastatic mechanisms under evolutionary pressure.

Sampling

The MHCC97L and the HCCLM6 Cells were harvested at logarithmic growth phase. DNA was extracted using DNA Extraction Kit (Agilent) for DNA sequencing and chip analysis. RNA was extracted using TaKaRa MiniBEST Universal RNA Extraction Kit for transcriptome analysis.

Whole-genome and whole-exome sequencing

The library preparation was performed as previously described. Briefly, paired-end sequencing libraries were prepared on 50 ng of cell genomic DNA using Illumina DNA Library Preparation Kit for whole-genome sequencing, and single-end sequencing libraries were prepared on 50 ng of captured DNA using Agilent SureSelect Human All Exon V5 Kit for whole-exome sequencing. The sequencing was performed on Illumina HiSeq 2000 instruments for WGS and SOLiD for WES. For the WES, the average depths of coverage were 88% and 84.5%, respectively for MHCC97L and HCCLM6, with 87.76% and 80.95% of their coverage larger than 30 X, and 54.66% and 50.28% of their coverage above 50 X (Supplementary Table 10).

Mutation calling

Raw sequencing reads of WES were aligned against human reference build GRCh37 using bwa (v0.7.11) [15]. The duplicates were removed by Picard (v1.4.5). And the base quality was recalibrated using the Genome Analysis Toolkit (GATK v4.) [16]. The two cell lines were derived from a liver cancer patient whose normal tissue was not sampled then, leaving us no control in the workflow. We therefore used two approaches for mutation calling. The first approach was to use canonical bcftools (v1.1.0) function to call mutations on the SAMtools (v0.1.6) [17] piled files in the reference of genome GRCh37, at a minimum coverage of 70X, mutation inference quality larger than 30 and mutated reads count cut off at three. Then we discarded the unique mutations for each cell line and retained the 980 of their overlapping mutations. Except for being entailed within the Catalogue of Somatic Mutations in Cancer (COSMIC) database (v70) [18], the SNVs were further filtered out if found in any of 1,000 Genomes Project (v3), duplicate gene flanking regions (v68), ENCODE DAC and Duke Mappability Consensus Excludable Regions (High Signal blacklist comprising poorly mapping reads and repeat regions), or the Fuentes’ suggestion list of likely false positive variants for exome sequencing [19]. Afterwards, the mutations with allele frequencies being higher than 50 percent in both cell lines were discarded, as it would be arbitrary to consider if the allele is wild type or mutant. Ultimately, 377 potential common mutations of the two cell lines were obtained. In the second approach we employed Mutect (v2.0) to predict post metastasis somatic mutations in each cell line with the other being the control. For specific somatic mutations in each cell line, VAF large than 0.1 was applied to filter the SNVs. One superior feature of Mutect is that it sensitively recognizes the situation where the tumor sample is believed to have mutated reads as well as the control sample, though with a much lower number. In usual circumstances, such mutations (labeled as Alterations in Normal or Germline Risk in the filters) are discarded as they cast an uncertainty to the somatic identity of the mutations. We, however, kept these mutations and regarded them as putative expanded mutations in the course of metastatic event. As concerned for expansion purposes, we determined that the Log 2 transformed variant allelic frequency ratio should outnumber 1, in addition to the filters that coverage larger than 30 and lower variant allelic frequency larger than 5 percent. 70 mutations were obtained from calling mutations using the HCCLM6 as the tumor and only 1 mutation was left when doing it on the contrary way by using the MHCC97L as the tumor. The SNVs were functionally annotated by ANNOVAR [20]. Nonsynonymous, stop-loss, stop-gain and splice-site SNVs (based on RefGene annotations) were considered to be functional. And Polyphen was used to predict functional influences of them. A polyphen score above 0.7 was defined as damaging.

Copy number profiling

To investigate the copy number variations, we performed Affymetrix SNP6.0 array-based cytogenetics copy number assay on the cells. According to the Affymetrix provided protocol, 500 ng genomic DNA from each cell line was prepared and hybridized to the arrays. Staining and washing were carried out on a Fluidics Station 450 (Affymetrix) and scanning with a GeneChip scanner 3000 7G system (Affymetrix). CEL files were generated by Affymetrix Command Console version 3.0 and imported into Genotyping Console (v 4.1 Affymetrix) to perform copy number and LOH analysis using regional, GC correction and the CN baseline file generated from 270 HapMap samples. The resulting CNCHP-files were used for genome position segmentation by R package DNAcopy from Bioconductor.

Next, we cross validated the copy number patterns of the two cell lines via applying the suite of HMMcopy (v1.14.0) on the WGS data. Briefly, the coverages were initially corrected with the GC and mappability bias of the reference genome. Then the corrected signals were segmented using Hidden Markov Model to yield an estimate of the copy number events. 7675 and 10526 autosomal segments were inferred respectively in the MHCC97L and the HCCLM6. After parameter optimization, the segment counts were reduced to 1669 and 1819.

Genomic Rearrangements

Mapped BAM files from WGS were used to call structural variants with Delly (v0.5.5) [21] at a minimum median mapping quality of 20 and a paired-end supporting reads cut-off of five. Then we used the variant reads and normal reads flanking the putative breakages to infer the frequencies of the structural mutating events. We regarded the SVs as the same event if both their corresponding breakage sites is within 5000 bp to each other’s flanking region. To identify genes affected by the genomic rearrangements, bed files were generated for each sample from deleted regions, and 1 kb region flanking the breakpoints from inversions, inter-chromosomal translocations, and tandem duplications. The resultant bed files were examined with Homer annotation. Structural events were visualized using sv tools and R package Circos (v0.67-4) [22].

Methylation microarray and analysis

Genomic global methylation was assessed on Illumina Human Methylation 450k Chip using 500 ng of input genomic DNA of each cell line. All digital analysis was performed in R statistical environment (v3.2.1). Intensities were converted to beta-values from the IDAT files using the wateRmelon package (v1.6.0) from Bioconductor (v3.0). Quality control was performed with the assistance of minfi package (v1.12.0). By default, raw data are filtered for probes with a detection P > 0.01 in at least one sample.

mRNA sequencing and analysis.

RNA from the two cell lines were extracted and purified with XX Kit. Libraries were made with these RNAs using mRNA-seq Library Prep Kit for Illumina. Then 75 base single-end sequencing was performed on the libraries on an Illumina HiSeq 2000 platform. Sequence reads were aligned to human genome version hg19 using TopHat 2.0. Normalized expression profiles with FPKM values were generated based on the accepted BAM files using the cufflinks suite. For gene expression comparisons, the human liver tissue RNA expressions were obtained from the HPA human RNA-Seq project via NCBI gene database. The HepG2 expression data was obtained from GEO dataset via accession number GSE88089

Subclonal reconstruction

R package Sciclone was used to characterize coexisting subpopulations (SPs) in the cell lines based on the point mutations and CNs with default settings. R package pracma was used to reconstruct the phylogenetic tree.

The mutational landscape of both MHCC97L and HCCLM6 cell lines

To explore the relationship between MHCC97L and the HCCLM6 cell populations, we first examined their karyotypes to uncover their chromosomal abnormality. The conventional karyotype analysis revealed that MHCC97L showed the heterogeneity with various forms of abnormal morphologies of chromosomes 1, 4, 6 and16, which exhibits at least three types of karyotypes (Supplementary 1), whereas HCCLM6 cells are relatively homogenous in the karyotypic appearances. This phenomenon that HCCLM6 cell populations are less diversified than the MHCC97L cells suggested that HCCLM6 cells could originate from a given cell subclone within MHCC97L cell populations. In addition, multiplex-fluorescence in situ hybridization (M-FISH) analysis indicated that the chromosomal karyotypes of the MHCC97L and HCCLM6 cells exhibited the dramatic aneuploidy and morphology abnormities (Supplementary Figure 1). Generally, the karyotypes of both cell lines are composed of 53 and 55 chromosomes respectively, with similar chromosomal aberrations. Chromosomal translocations between chromosomes 1 and 21, 4 and 8, 9 and 7, 9 and 2, 14 and 22 were common in the two cell lines, whereas HCCLM6 cells showed extra translocations between 1 and 8, 16 and 3, as well as more numbers of chromosomes 7 and 16. This implied that HCCLM6 cells could undergo the featured genetic events during metastatic clonal evolution.

To test the hypothesis, we further performed whole genomic analysis on the two cell lines, including whole genome sequencing (WGS), whole exome sequencing (WES) and Affymetrix single-nucleotide polymorphism (SNP) genotyping as well as methylome profiling with Illumina Infinium HumanMethylation450 BeadChip and the transcriptome profiling with RNA sequencing (RNA-Seq), according to the workflow illustrated in Supplementary Figure 2. And then genomic aberrations including single nucleotide variant (SNV), copy number variation (CNV) and structural variant (SV) based on WGS, WES and SNP genotyping were called in both cell lines (see Methods).

Briefly, through SNV calling and filtering on the WES data, we obtained 335 overlapping mutations, as well as 165 and 160 mutations exclusively in MHCC97L and the HCCLM6 cells, respectively (Supplementary Table 1), of which 347 were nonsilent mutations. CNV data was obtained based on Affymetrix SNP 6.0 genotyping (Supplementary Table 2). Besides, based on the WGS data, whole genomic structural aberrations, including CNV (Supplementary Table 3) and genomic rearrangements (Supplementary Table 4), were obtained. The comparison of the two above CNV data showed that Pearson’s correlation coefficient is 0.91 (P < 2.2e-16) between them, indicating the high data reliability. Meanwhile, to supplement the information about the epigenetic status and transcriptional activities of the two cell lines, DNA methylation profiles and the gene expression profiles were yielded (Supplementary Tables 5 and 6). The genome-wide representation of these results depicted with Circos plots (Figure 1), in general, revealed that the two cell lines presented similar genomic aberrations and highly homologous characteristics, supporting that they share common ancestry. However, some distinctive features might drive the metastatic clonal evolution within HCCLM6 cell populations.

A minor subclone in the MHCC97L cells became dominant in the HCCLM6 cells

To study the heterogeneities and phylogenetic relationships in MHCC97L and HCCLM6, here we performed clonal analysis on the two cell lines based on the SNVs and CNVs (see Methods). The mutations were deconvoluted into 3 subclones (denoted as 1, 2 and 3) within the MHCC97L cells (Figure 2A), and 2 subclones (denoted as A, B) in the HCCLM6 cells (Figure 2A). Then we counted the mutations in each subclone. Interestingly, the minor clones, subclone 3 in the MHCC97L and subclone B in the HCCLM6, are mainly characterized with their specific mutations, whereas the shared mutations were expanded in subclone 1, 2 and 3 in the MHCC97L, and majority in subclone A in the HCCLM6 (Figure 2B). The subclonal structure of the cell lines complied with notion that the MHCC97L could be more heterogenous. In addition, the results indicated that the subclone A in the HCCLM6 could be originated from subclone 3 in the MHCC97L. To confirm the implication, next, we performed phylogenetic reconstruction based on the mutated alleles. It was indicated that the HCCLM6 are derived from the same node with subclone 3 in the MHCC97L (Figure 2C), which consolidate the notion that the ancestry minor subclone 3 in MHCC97L, when it first emerged, with all the MHCC97L and HCCLM6 shared mutations, and without the current MHCC97L unique mutations, is the ancestor of the subclone A in the HCCLM6 (Figure 2D). It was noticed that c.2156C>T in DMBT1 was implicated as early mutation in the evolutionary history as it was in subclone 1 and A respectively in the MHCC97L and HCCLM6 (Supplementary Table 1). The newly formed functional damaging nonsilent mutations in the MHCC97L includes those influencing DDX59, SETD5, CXCR6, CDK14, etc. (Supplementary Table 1). Meanwhile, newly formed functional damaging nonsilent mutations in the HCCLM6 includes those influencing EFHC2, BCORL1 etc. (Supplementary Table 1). The nonsilent mutations that were specifically expanded from the subclone 3 in the MHCC97L to the subclone A in the HCCLM6 included c.1061T>G in SMAD5 and c.1373A>G in RNF169 (Supplementary Table 1 and Supplementary Figure 3), both of which were predicted as functional damaging.

MHCC97L and HCCLM6 cells concurrently harbored amplification of CAV1, CAV2 and MET

Due to the HCCLM6 cells originating from the subclone of MHCC97L, the common genetic events shared by the two cell lines could be initial drivers responsible for metastasis.

Next, the unstable genome integrities of the two cell lines were evaluated by SV and CNV analyses. For SV, we observed 24 concordant inversions and 30 concordant translocations in the two cell lines. Moreover, the MHCC97L harbored 44 inversions and 15 translocations exclusively, whilst the HCCLM6 had 8 unique translocations and none unique inversions (Figure 1). Annotation of the SV events revealed that they did not flank gene coding regions.

The CNV analysis using SNP6.0 array revealed that the average autosomal ploidy for the MHCC97L and the HCCLM6 were respectively 2.03 and 2.06. The MHCC97L includes 121 CNV gain (average segment size 3.2 Mb) and 121 CNV loss (average segment size 1.6 Mb); the HCCLM6 includes 128 CNV gain (average segment size 3.1 Mb) and 145 CNV loss (average segment size 1.6 Mb). As shown in the global track of the genomic CNVs, the MHCC97L and HCCLM6 have undergone analogous copy number altering events (Figure 3A), with typical amplification on chromosome 7 (Figure 3B) and the deepest deletions on chromosome 8 and 2 (Figure 3C and 3D). In consistence with previous studies, we have observed chromosome 8p loss in the two cell lines, especially the region harboring the genes of LOC392196, USP17L7, ZNF705D, FAM66D and USP17L2 (Figure 3C). Chromosome 8p loss is widely studied and testified to associate with tumorigenesis and metastasis as it harbored several tumor suppressing genes. Interestingly, the 8P loss is slightly deeper in the HCCLM6 than in the MHCC97L, which is consistent, as well as most of the genome, in the WGS and Chip revealed results (Supplementary Figure 4A). In despite of the concordance, the both ends of chromosome 2, as well as the head of chromosome 3p, chromosome 15, are differed in the two cell lines in their CN states.

The genes with copy number larger than 6 or lower than 1 was listed in Supplementary Table 7, among which The chromosome 7q31.2 region flanking CAV1, CAV2 and MET was the mostly amplified region in the whole genome and therefore constitute the prominent peaks in the CNV profiles in the two cell lines (Figure 3B). MET amplification is a significant recurrent CNV event in the TCGA HCCs according to cBioportal. The proto-oncogene, have amplified along with CAV1 and CAV2 to reach more than 10-13 copies in the MHCC97L and HCCLM6. In contrast, the most deeply deleted gene in the two cell lines lied in the segment of Chr8:11968639-12017716, followed by Chr2:191833566-192049769, where tumor suppressor gene STAT1 is located (Figure 3D). We compared the gene expressions of these most amplified and deleted genes in the two cell lines, as well as in the normal liver tissue and the HepG2 cell line, which do not report such alterations. As shown in Supplementary Figure 4B, the dose effect of the gene copies seemed obvious. In contrary to the corresponding expressions in the normal liver and HepG2, CAV1, CAV2, MET were largely overexpressed in the two cell lines, especially in the HCCLM6, while STAT1 was utterly silenced (Supplementary Table 8). Amplification and overexpression of MET are reported to be associated with cell survival and invasiveness in multiple human cancers, especially HCC. Therefore, it is highly likely that MET amplification was a driver mutation that propelled the malignancy and metastasis in the early tumor phase of the MHCC97.

EMT associating genes were remarkably elevated in HCCLM6

To further investigate the dose effects of CNV events on gene expressions, we performed intra- and inter-cell comparisons. Firstly, for the intra- cell comparison of gene expressions, we assigned the genes into four groups according to their copy numbers and then observed a significant overall increase (P= 1.3e-5) of gene expressions along with the increment of gene copies in both cell lines (Supplementary Figure 4C), suggesting that increasing gene dose had generally elevated the transcription activities.

Next, for inter-cell comparisons, we compared the expressions of the genes that have different copy numbers in the MHCC97L and the HCCLM6. Interestingly, the group of genes that had more DNA copy numbers in the HCCLM6 had up-regulated expressions compared to their correspondents in the MHCC97L (Supplementary Figure 4D). In contrast, the genes that had concordant copy numbers were highly expressed in the MHCC97L (P=2.7e-4) and the genes that had less copy numbers in the HCCLM6 had no significant differences in their total expression distributions in both cell lines. The result showed that, in despite of dose effect of copy numbers, the global gene transcriptions in the HCCLM6 are more active.

Among the genes that had differed copy numbers between the two cell lines, the top ranked genes that had gained copy numbers and upregulated transcriptions in the HCCLM6 in comparison to the MHCC97L are CAV1, CAV2 and MET. Although these genes were amplified in both cell lines, the HCCLM6 had 2 more copies and had an elevation of more than 100-fold in the expressions. Previous studies reported that the Caveolins had significantly facilitate portal vein invasion of HCC. They can trigger the EMT process by up-regulating Vimentin and inhibit E-cadherin. Also, it increased secreted MMPs. Therefore, we examined the expression of these genes in the cell lines. It was found that N-cadherin, which helped the tumor cells to interact with stromal cells, is 23-fold higher expressed in the HCCLM6. In contrast, E-cadherin that guarantees tight junction is lowly expressed in both cell lines. Moreover, the secretory protease of MMP7 is 6-fold higher expressed in the HCCLM6. There is also implication that CAV1 and CAV2 favors metastasis via exome secretion pathway as they were in the exosome cargoes of metastatic HCCs.

In addition, in the substantial genes that are slightly differed between the two cell lines, we noticed that EMT transcription factor TWIST2 was among it; it had a shallow deletion in the MHCC97L but not in the HCCLM6. Then we compared the expressions of these three genes of the cell lines with the liver tissue and the HepG2. The trend of the expressions of the genes reflect to their malignancy states (Supplementary Table 8). Their expressions were elevated in the MHCC cells than in the liver and the HepG2. Besides, they were largely up-regulated in the HCCLM6 than in the MHCC9L.

DLEC1 is aberrantly methylated in the two cell lines

To investigate the DNA methylome of the two cell lines, we performed HK450 Array hybridization analysis. As illustrated in Supplementary Figure 5 and Supplementary Table 4, the probe intensities had similar distribution in the two cell lines. In addition, we observed that larger than 99 per cent of the probes had approximate intensities (differences within 0.1), which indicated that the two cells are similar in the landscapes of their genomic methylations. The two cell lines had concurrent high densities of methylation in the promoter region of tumor suppressor gene DLEC1. As shown in Figure 4A, the promoter region of DLEC1 is occupied with a CpG island and its N-terminus shore region, which are densely marked with high methylations intensities. Besides, the methylations in the HCCLM6 are slightly higher than in the MHCC97L. The DLEC1 expressions is down-regulated by 19- and 36- fold change in MHCC97L and HCCLM6 respectively than in the normal liver tissue.

In the purpose of screening differentially methylated genes, we chose the genes that had at least deviated methylated CpG sites that exhibited more than three tandem sites that had intensity differences larger than 0.3 in the promoter region (Supplementary Table 9). We selected the most heavily marked genes and aligned them with elevated expressions, among which top ranked genes are TP6V0E2, ZC3HAV1L, SGCE, NIPA1, YTHDC1, HUWE1, GRB10 and SNRPN that are relatively hypomethylated in the HCCLM6.Several encompassed genes are strong biological candidates for tumorigenesis and metastasis. Specifically, YTHDC1, HUWE1, VIM and GRB10, which are elevated with 27.5-, 8- and 7.6- fold change in the HCCLM6, are involved in mRNA splicing, tumor suppressor ubiquitination and tyrosine signaling respectively (Figure 4B).

Next, we examined the pathways of the cell lines via performing enrichment analysis on the transcription profiles. The result showed that the MHCC97L exceled in the chromatin package, cell cycle and proliferation pathways (Supplementary Figure 6). In contrast, the HCCLM6 transcription activities are more enriched in the protein folding, exome secretion processes, which is explicitly in consistence with the implications of the up-regulated expressions of CAV1 and CAV2.

HCCLM6 exceled via reshaping the clonal structure rather than obtaining novel alterations

It is worth mentioning that the most implicative metastatic favoring alterations are presented concordantly in the MHCC97L and HCCLM6. For instance, in addition to the well mentioned MET, CAV1, CAV2 amplification, the deletion of STAT1 implicates favors to metastasis as its high expression is reported to suppresses angiogenesis in human glioma cells. Moreover, the hypermethylated DLEC1, as a tumor suppressor, is reported to be deleted or aberrantly methylated in multiple cancer forms. All the superiority of HCCLM6 in metastasis are likely to be based on these same materials, merely with minor elevations in the damaging extents. Specifically, the further amplification of MET,CAV1, CAV2 in the HCCLM6 had resulted in larger gene doses and remarkable transcriptionally elevations. Additionally, the DLEC1 promoter methylation in the HCCLM6 is slightly lower and so is the expression. Moreover, the promoter region of HUWE1 is mediately methylated in the MHCC97L and become largely erased in the HCCLM6. HUWE1 encodes a protein that functions as an E3 ubiquitin ligase that helped degradation of the TP53 tumor suppressor, core histones, and DNA polymerase beta. Given the fact that clonal analysis revealed that the HCCLM6 is likely to be formed via expansion of an early minor subclone like that in the MHCC97L, these mentioned changes in the extents of damaging alterations suggested that the HCCLM6 may become aggressively metastatic via reshaping the clonal structure to push its evolution in the course towards the metastasis favoring path. To conclude, we proposed a schema chart to illustrate the relationships between the two cell lines with the aforementioned important mutations under the clonal expansion proposition (Supplementary Figure 7).

The complexities of the timing, sites and phenotypic shifts of metastases are persistently mysterious due to its insidious camouflage and evolving paradigm. The cancerous cells exploit multiple distinct biological mechanisms to advance through their courses of metastasis cascades. Some testified events encompass the traditional theory of epithelial-mesenchymal transitions (EMT), extracellular matrix remodeling, tumor-stroma interactions and the newly proposed pioneering exosome secretion, as well as the competent state initialized via chromatin opening [2-4, 14]. Although the cell lines of MHCC were well established for metastatic analyses, the driving genomic alterations behind the scene are largely unknown.

It is accepted that transformed cells, with undermined integrities, undergo proliferative progression and successively acquire novel aberrations to generate genetically diverse subclonal architectures. Driver mutations contributing to the cell lines’ differentiating metastatic performances should have been carried on by their correspondent progenitor cells, which means existence preceded the time course of the expanding. The result in this study suggested a scenario that the emergence of certain mutations in the MHCC97 with low frequencies, which was stable inherited to the MHCC97L, could have endowed the cells a fair low metastatic rate as attested in the previously reported orthotopic implantation experiments. Subsequently, after artificial selection through the recycling processes of implantation, metastasis and dissection, the offspring of the metastatic minorities in the MHCC97 may accumulate high frequencies of these early metastasis driving mutations and therefore become dominant in population, making the entity a highly metastatic one.

Based on the results, we propose that one strategy whereby HCCLM6 acquired the metastatic capabilities might be its enhancement of EMT associating activities. An interesting contradiction in the cell lines is that although the HCCLM6 is more invasive, it is inhibited comparing to MHCC97L in the perspective of cell cycle and proliferation, as suggested by the methylation regulated genes and the RNA-Seq revealed global transcription profiles. However, cell division is limited in the mesenchymal state than in the epithelial state. And the expression of N-Cadherin and MMP7 in the HCCLM6 indicated the onset of epithelial to mesenchymal transitions. There are also implications that the HCCLM6 had EMT favoring events, such as its higher gene dose of TWIST2 and the corresponding transcriptional upregulations. TWIST2 is an EMT orchestrating transcription factor that would enhance the cells’ capability to associating with stroma cells during invasion process. Therefore, the phenomenon makes sense considering the foreseeable local tradeoff between cell growth and invasion. Another strategy that worth mentioning is the CAV1 and CAV2 associated pathways. In addition to its associations with EMT, CAV1 and CAV2 was also reported to play metastasis facilitating roles via extracellular vesicles transport and exosome secretion according to previous studies, which were identified to be enriched in the HCCLM6 via the RNA-Seq data.

To conclude, considering the main causes to explain the metastasis discrepancy between the MHCC97L and HCCLM6 include, this study revealed that, 1) MET, CAV1 and CAV2 amplifications and STAT1 deep deletions are likely to be the main drivers for the malignancy and low metastatic capability for MHCC97L, 1) a minor subclone is expanded to result in the highly invasiveness in the HCCLM6 and dysfunctional DMBT mutation is expanded, and 3) aberrant methylation and CNVs in the HCCLM6 resulted in the high expressions of genes that favor the EMT and exome secretion in tradeoffs with cell growth to presumably accelerate metastasis.

HCC: Hepatocellular carcinoma; MHCC: human Metastatic Hepatocellular Carcinoma; HBV: hepatitis B virus; HCV: hepatitis C virus; WGS: whole genome sequencing; WES: whole exome sequencing; SNP: single-nucleotide polymorphism; SNV: single nucleotide variant; CNV: copy number variation; SV: structural variant; M-FISH: multiplex-fluorescence in situ hybridization; EMT: epithelial-mesenchymal transitions.

Ethics approval and consent to participate

Not applicable

Consent for publication

Not applicable

Availability of data and material

The mouse next-generation sequencing data used in the manuscript can be downloaded from the database of NCBI under accession number of SUB6534381

Competing interests

The authors declare that they have no competing interests.

Funding

This work was supported by the China National Science and Technology Major Project for Prevention and Treatment of Infectious Diseases grant (2017ZX10203207, to Z.G.H.)

Authors' contributions

Q L, YS and NW drafted the manuscript. ZN L and ZG H edited the manuscript. All authors read and approved the final manuscript.

Acknowledgements

Not applicable

El-Serag HB, Rudolph KL: Hepatocellular carcinoma: epidemiology and molecular carcinogenesis. Gastroenterology 2007, 132(7):2557-2576.
Fraser M, Sabelnykova VY, Yamaguchi TN, Heisler LE, Livingstone J, Huang V, Shiah YJ, Yousif F, Lin XH, Masella AP et al: Genomic hallmarks of localized, non-indolent prostate cancer. Nature 2017, 541(7637):359-+.
Muller A, Homey B, Soto H, Ge NF, Catron D, Buchanan ME, McClanahan T, Murphy E, Yuan W, Wagner SN et al: Involvement of chemokine receptors in breast cancer metastasis. Nature 2001, 410(6824):50-56.
He M, Qin H, Poon TCW, Sze SC, Ding XF, Co NN, Ngai SM, Chan TF, Wong N: Hepatocellular carcinoma-derived exosomes promote motility of immortalized hepatocyte through transfer of oncogenic proteins and RNAs. Carcinogenesis 2015, 36(9):1008-1018.
Heldin CH, Landstrom M, Moustakas A: Mechanism of TGF-beta signaling to growth arrest, apoptosis, and epithelial-mesenchymal transition. Curr Opin Cell Biol 2009, 21(2):166-176.
Hugo H, Ackland ML, Blick T, Lawrence MG, Clements JA, Williams ED, Thompson EW: Epithelial--mesenchymal and mesenchymal--epithelial transitions in carcinoma progression. J Cell Physiol 2007, 213(2):374-383.
Pollard JW: Tumour-educated macrophages promote tumour progression and metastasis. Nat Rev Cancer 2004, 4(1):71-78.
Denny SK, Yang D, Chuang CH, Brady JJ, Lim JS, Gruner BM, Chiou SH, Schep AN, Baral J, Hamard C et al: Nfib Promotes Metastasis through a Widespread Increase in Chromatin Accessibility. Cell 2016, 166(2):328-342.
Sun FX, Tang ZY, Lui KD, Ye SL, Xue Q, Gao DM, Ma ZC: Establishment of a metastatic model of human hepatocellular carcinoma in nude mice via orthotopic implantation of histologically intact tissues. Int J Cancer 1996, 66(2):239-243.
Tian J, Tang ZY, Ye SL, Liu YK, Lin ZY, Chen J, Xue Q: New human hepatocellular carcinoma (HCC) cell line with highly metastatic potential (MHCC97) and its expressions of the factors associated with metastasis. Br J Cancer 1999, 81(5):814-821.
Li Y, Tang ZY, Ye SL, Liu YK, Chen J, Xue Q, Gao DM, Bao WH: Establishment of cell clones with different metastatic potential from the metastatic hepatocellular carcinoma cell line MHCC97. World J Gastroenterol 2001, 7(5):630-636.
Li Y, Tang ZY, Ye SL, Liu BB, Liu YK, Chen J, Xue Q: Establishment of a hepatocellular carcinoma cell line with unique metastatic characteristics through in vivo selection and screening for metastasis-related genes through cDNA microarray. J Cancer Res Clin 2003, 129(1):43-51.
Li Y, Tian B, Yang J, Zhao L, Wu X, Ye SL, Liu YK, Tang ZY: Stepwise metastatic human hepatocellular carcinoma cell model system with multiple metastatic potentials established through consecutive in vivo selection and studies on metastatic characteristics. J Cancer Res Clin 2004, 130(8):460-468.
Baca SC, Prandi D, Lawrence MS, Mosquera JM, Romanel A, Drier Y, Park K, Kitabayashi N, MacDonald TY, Ghandi M et al: Punctuated Evolution of Prostate Cancer Genomes. Cell 2013, 153(3):666-677.
Li H, Durbin R: Fast and accurate long-read alignment with Burrows-Wheeler transform. Bioinformatics 2010, 26(5):589-595.
McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, Garimella K, Altshuler D, Gabriel S, Daly M et al: The Genome Analysis Toolkit: a MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res 2010, 20(9):1297-1303.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, Marth G, Abecasis G, Durbin R: The Sequence Alignment/Map format and SAMtools. Bioinformatics 2009, 25(16):2078-2079.
Forbes SA, Beare D, Gunasekaran P, Leung K, Bindal N, Boutselakis H, Ding MJ, Bamford S, Cole C, Ward S et al: COSMIC: exploring the world's knowledge of somatic mutations in human cancer. Nucleic Acids Res 2015, 43(D1):D805-D811.
Fajardo KVF, Adams D, Mason CE, Sincan M, Tifft C, Toro C, Boerkoel CF, Gahl W, Markello T, Sequencing NC et al: Detecting false-positive signals in exome sequencing. Hum Mutat 2012, 33(4):609-613.
Wang K, Li MY, Hakonarson H: ANNOVAR: functional annotation of genetic variants from high-throughput sequencing data. Nucleic Acids Res 2010, 38(16).
Rausch T, Zichner T, Schlattl A, Stutz AM, Benes V, Korbel JO: DELLY: structural variant discovery by integrated paired-end and split-read analysis. Bioinformatics 2012, 28(18):I333-I339.
Krzywinski M, Schein J, Birol I, Connors J, Gascoyne R, Horsman D, Jones SJ, Marra MA: Circos: An information aesthetic for comparative genomics. Genome Research 2009, 19(9):1639-1645.

Download PDF

Version 1

posted

You are reading this latest preprint version

Genomic features reveal the clonal evolution in human hepatocellular carcinoma cell lines MHCC97 with stepwise metastatic potentials

Status:

Version 1

Abstract

Figures

Background

Methods

Results

Discussion

Conclusions

Abbreviations

Declarations

References

Supplementary Files

Status:

Version 1