This article is mentioned in:

Introduction

Hepatocellular carcinoma (HCC) is one of the most lethal malignancies worldwide due to its aggressive properties and limited therapeutic options. Thus, one important aim of the present study was to identify biomarkers for predicting therapeutic response and prognosis, which may provide a deeper understanding of the tumor biology (1,2). To date, numerous molecular events, including those at the epigenetic, pre-transcriptional, transcriptional, translational and post-translational levels, have been identified and validated as important prognostic biomarkers (3), such as the Ki-67 protein and TP53 gene mutations. However, the functional significance of cancer-specific alternative splicing (AS) events is largely unexplored. AS of pre-mRNA enables a gene to produce multiple protein isoforms, usually with distinct functions. This mechanism generates transcriptomic and proteomic diversity in different tissues and cell types (4,5), however, in some disease states, for example in cancerous conditions, the damaged or diseased cells take this advantage to produce aberrant protein isoforms that contribute to tumorigenesis (6).

In the last decade, growing data has indicated that the dysregulation of AS is a crucial event in carcinogenesis and tumor progression (3,7–9). Genes with aberrant AS events are involved in almost every aspect of cancer, including proliferation, differentiation, cell cycle control, metabolism, apoptosis, motility, invasion and angiogenesis (10). Although several studies have revealed multiple tumor-associated AS variants in a wide variety of cancers (10–14), there is still much research to be conducted at the splicing level in HCC. In HCC, various aberrant AS events were found to promote the generation of oncogenic isoforms, whereas tumor suppressors were self-inactivated by aberrant AS (15). As a result, the identification and validation of isoforms as novel biomarkers at the splicing level in HCC may be of great significance for revealing novel drug targets for therapeutic intervention.

In the present study, we analyzed the published RNA-Seq datasets of three paired HCC and adjacent non-tumor tissues and summarized the aberrant gene expression and splicing patterns in HCC samples. We focused on the putative novel biomarkers of these three aberrant AS events in HCC, and their association with clinical features and prognosis was further investigated.

Materials and methods

RNA-Seq data download and processing

The RNA-Seq data of three paired HCC patients were downloaded from NCBI GEO database with the accession number of GSE33294 (16). Read alignment was performed using TopHat2 (17) with the known reference genes from the GENCODE v22 annotation (http://www.gencodegenes.org/human/release_22.html). Gene expression was analyzed using Cufflinks v2.2.1 (http://cole-trapnell-lab.github.io/cufflinks/) and significantly differential genes were identified by the Cuffdiff module at the FDR of 5% (18). Alternative splicing analysis was carried out using the MATS software v3.0.8 (http://rnaseq-mats.sourceforge.net/) with the parameters of ‘-t single-len 58 -a 8 -c 0.0001 -analysis U’. Only the significantly differential events of the SE type were used for subsequent analysis.

Clinical samples and ethical statement

Forty-five paired tumor and normal liver tissues from HCC patients, with an average age of 45 years old, including 28 males and 17 females, were collected from April 2010 to June 2014 at the China-Japan Union Hospital, Jilin University. All patients consented to the use of their clinical specimens for medical and biological study. The present study was approved by the Clinical Research Ethics Committee of China-Japan Union Hospital, Jilin University (2015-wjw002).

RT-qPCR validation and inclusion level (IL) calculation

Fresh samples were frozen in liquid nitrogen and total RNA was extracted using a TRIzol Isolation kit. cDNA was synthesized using the PrimeScript™ RT Master Mix (Takara Bio, Inc., Otsu, Japan). RT-PCR was performed on the Illumina ECO system (Illumina, Inc., San Diego, CA, USA). ILs of skipped exons in every clinical sample were estimated by the ΔΔCq approach (19). Each splicing event used two paired primers, one representing the inclusion isoform expression and the other representing the total expression. The inclusion ratio was calculated using the following equation: IL=2-(Inclusion ΔCT−TotalΔCT). Thus, the ‘InclusionΔCT-TotalΔCT’ was used to denote the IL of each skipped exon.

Association analysis between AS and clinical features

Clinical features for the association analysis included age, sex, viral status, cell differentiation, tumor size, vascular invasion, alpha-fetoprotein (AFP) level, metastasis and recurrence, all of which were converted to the quantitative data. The correlation between these quantitative clinical features and splicing changes was calculated using the maximal information coefficient (MIC) in the MINE toolkit (20) with the parameter of ‘-allPairs’. The confidence for various MIC values was evaluated by the P-value table at the sample size of 45.

Statistical analysis and visualization

Heatmaps of expression and splicing were plotted by the ‘pheatmap’ package in R language with the parameter of ‘scale’=‘row’. Three DS genes were visualized using the sashimi plotting function in IGV software (21). Student's t-test was used in the differential analysis of RT-PCR values using GraphPad Prism version 5.0 (GraphPad Software, Inc., La Jolla, CA, USA). Fisher's exact test was applied to the association test of contingency table analysis. The network of AS and clinical features was visualized using Cytoscape software v3.3.0 (http://cytoscape.org). Enrichment analysis was performed using the ‘ClusterProfiler’ package (22) in R language.

Results

Comparison between gene expression and splicing profiles of HCC

We reanalyzed the RNA-Seq datasets from three previously published paired HCC tumors and adjacent non-tumor tissues (16). After reads alignment and expression quantification, a total of 13,327 protein-coding genes were identified with FPKM≥1. Of these, we detected 1,246 upregulated genes and 1,197 downregulated genes in the HCC group at a false discovery rate (FDR) of 5% using the Cufflinks software (18). The expression pattern for these differentially expressed (DE) genes in three paired samples is presented in Fig. 1A, and it was revealed that the non-tumor and tumor samples were clearly separated and the size of upregulated and downregulated gene clusters were nearly equal. Additionally, there were similar expression patterns between three samples in the same group of tumor or non-tumor samples.

Figure 1. Expression and splicing profiles of HCC tumors and non-tumors. (A) Unsupervised clustering of differentially expressed genes between HCC tumors and non-tumors. The FPKM values of genes were inputted. T448, T473 and T510 were the HCC tumor sample numbers. N448, N473 and N510 were the non-tumor sample numbers. (B) Unsupervised clustering of differentially spliced events between HCC tumors and non-tumors. The PSI values of splicing events were inputted. (C) Overlap between differentially expressed genes and differentially spliced genes. HCC, hepatocellular carcinoma; FPKM, fragments per kilobase of transcript per million mapped reads; PSI, percent-spliced-in; DE, differentially expressed; DS, differentially spliced.

Besides gene expression changes, AS, as an important post-transcriptional modification, also caused differences between tumor and non-tumor tissues. The splicing differences were analyzed using MATS (23,24) for five common splicing types, including skipped exon (SE), retained intron (RI), mutually exclusive exons (MXEs), alternative 5′splice site (A5SS), and alternative 3′splice site (A3SS). The most frequent splicing type was the SE type, which had 197 splicing events with statistically significant differences. The percent-spliced-in (PSI) value, representing the fraction of the exon-inclusion variant, was used to estimate the splicing level. As in the expression patterns of DE genes, the hierarchical clustering for PSIs of differential skipped exons still illustrated the separation between the tumor and non-tumor samples (Fig. 1B). However, PSI-decreased exons were obviously more than PSI-increased exons, suggesting that HCC broke the splicing balance of exon inclusion and exclusion on the whole. The differentially spliced (DS) genes were fewer in number than the DE genes, and only 29% of the former overlapped with the latter (Fig. 1C). Moreover, the splicing specificity of samples also existed in the same group. These splicing-changed genes in HCC revealed a new layer of regulation in liver carcinogenesis.

Function and pathway enrichment of DE and DS genes

Next, an enrichment analysis was performed on DE and DS genes between HCC tumor and non-tumor samples to observe whether the splicing-induced changes were clustered in the main functions and processes of the liver and liver cancer. Top differentially-enriched terms of biological processes (P<0.05) between DE and DS genes are presented in Fig. 2A. It was revealed that the DE gene-specific-enriched terms included carboxylic acid catabolic process and cellular amino acid catabolic process. DS genes were primarily enriched in Gene Ontology (http://geneontology.org) terms involving acute-phase response, platelet degranulation and acute inflammatory response, which were identical to the core roles of the liver in immunological effects. Similarly, enrichment analysis of the KEGG pathway (https://www.genome.jp/kegg/pathway.html) revealed the enriched pathways of liver characteristics, such as cytochrome P450 metabolism, and complement and coagulation cascades in DS genes (Fig. 2B). These results indicated that the DS genes reflected the important functions and processes of the liver even though they were far fewer than the DE genes.

Figure 2. Enrichment analysis of expression and splicing profiles. (A) The enriched biological process terms of DE and DS genes. (B) The enriched KEGG pathway terms of DE and DS genes. The gradient color of the dots indicates the adjusted P-value of each enriched term. The size of the dots indicates the percentage of genes of each corresponding term in total differential genes. DE, differentially expressed; DS, differentially spliced.

Clinical sample validation for AS isoforms as novel candidate biomarkers

As aforementioned, the splicing changes were indispensable for investigating the differences between HCC tumors and non-tumors as well as exploring HCC biomarkers. Among the 197 DS events, we chose three AS events from both the PSI values and the functional importance of the related genes to validate their potential as novel splicing biomarkers in HCC. Reads alignments and coverages in three paired samples illustrated the splicing junctions of these three DS events (Fig. 3). We found that exon 7 of carcinoembryonic antigen-related cell adhesion molecule 1 (CEACAM1) and exon 2 of VPS29 revealed an increase in exon inclusion in HCC samples. CEACAM1 encodes carcinoembryonic antigen-related cell adhesion molecule 1, and the isoform with exon 7 inclusion is the longest transcript of CEACAM1 and plays a role as a co-inhibitory receptor in the immune response (25). VPS29 is a vacuolar protein sorting-associated protein and the exon 2-included isoform adds four amino acids to the protein initiation position. Moreover, exon 3 of ISOC2 revealed an increase in exon exclusion in HCC samples, the exclusion of which can alter the core isochorismatase domain of ISOC2.

Figure 3. Read coverage of three differentially spliced genes. (A-C) Sashimi plots revealed the sequencing read coverage for differentially spliced events of (A) CEACAM1, (B) VPS29 and (C) ISOC2 in all six samples. Sample numbers from top to bottom are N448, T448, N473, T473, N510 and T510. The number in each sashimi plot represents the read counts covering different exon junctions. CEACAM1, carcinoembryonic antigen-related cell adhesion molecule 1.

Although the splicing differences were identified by mRNA sequencing from only three paired samples, a large number of samples were employed to strictly validate the possibilities of splicing isoforms as biomarkers of HCC. The exon inclusion level (IL) differences were assessed by RT-qPCR in a total of 45 paired HCC tumor and adjacent non-tumor samples (Fig. 4). The ΔΔCq values of DS exons were significantly different between tumor and non-tumor tissues (P<0.01) and the splicing variation tendency was the same as the evidence from the RNA sequencing analysis. These splicing switches were expected to be closely associated with HCC pathology.

Figure 4. RT-PCR validation of three differentially spliced genes. RT-PCR of 45 paired HCC tumor and adjacent non-tumor samples revealed significant differences in IL for CEACAM1 exon 7, (A) VPS29 exon 2 (B) and ISOC2 exon 3 (C). N, non-tumors; T, tumors; HCC, hepatocellular carcinoma; CEACAM1, carcinoembryonic antigen-related cell adhesion molecule 1. ***P<0.0001. HCC, hepatocellular carcinoma; IL, inclusion level; CEACAM1, carcinoembryonic antigen-related cell adhesion molecule 1.

Association of DS with clinical features of HCC

Clinicopathological analysis of the three DS genes in 45 HCC patients was performed to evaluate which features were significantly associated with splicing switching. Contingency table analysis and Fisher's exact test revealed that the exon 7 IL of CEACAM1 was significantly associated with the viral status of clinical samples (Table I and Fig. 5A). HBV-/HCV-infected patients were more likely to express the isoform with exon 7 inclusion. The difference test also represented the significant difference of ΔIL, that is IL(tumor)-IL(non-tumor), between the HBV-/HCV- and non-infected patients (Fig. 5B).

Figure 5. Clinicopathological analysis of differentially spliced genes in 45 HCC patients. (A) Distribution of the IL changes of CEACAM1 exon 7 in HCC patients with distinct viral statuses. ‘IL(T)-IL(N)’ represents the IL change between tumor groups and non-tumor groups. (B) Difference test of the IL changes of CEACAM1 exon 7 between the viral-infected and non-viral patients. *P<0.05. (C) Association network of splicing or splicing changes with clinical features. The thickness of the edges indicates the MIC values of each paired node. ΔIL-, IL- and I- separately represent the IL change between HCC tumors and non-tumors, the IL in HCC tumors and inclusion isoform expression in tumors. -C, -I, -V separately represent CEACAM1, ISOC2 and VPS29. (D) Survival curves for CEACAM1 and its exon 7 inclusion isoform. (E) Scatter plot revealing the association between the survival time and ΔIL of CEACAM1 exon 7. (F) Survival curves for VPS29 and its exon 2 inclusion isoform. (G) Scatter plot revealing the association between cell differential stage and ΔIL of VPS29 exon 2. CEACAM1, carcinoembryonic antigen-related cell adhesion molecule 1.

Table I. Contingency table analysis of CEACAM1, VPS29 and ISOC2 splicing changes.

Table I.

Contingency table analysis of CEACAM1, VPS29 and ISOC2 splicing changes.

	CEACAM1				ISOC2			VPS29
Clinical features	Number	+	−	P-value	+	−	P-value	+	−	P-value
Age (years)
>50	25	23	2	0.6378	6	18	0.2696	22	2	1
≤50	20	16	3		2	17		18	2
Sex
Male	28	23	5	0.1414	4	24	0.4188	24	4	0.2797
Female	17	16	0		4	11		16	0
Viral status
Yes	36	34	2	0.0471	7	28	1	32	4	0.5687
No	9	6	3		1	8		9	0
Cell differentiation
Poorly	3	2	1	0.3188	0	3	0.8444	3	0	0.6662
Poorly-moderately	5	4	0		1	3		5	0
Moderately	23	21	2		4	18		20	2
Moderately-well	2	1	1		0	2		2	0
Well	4	3	1		1	3		3	1
Tumor size (cm)
>5	17	15	1	0.638	3	13	1	15	2	0.6337
≤5	28	24	4		5	22		25	2
Vascular invasion
Yes	9	8	1	1	2	6	0.6354	9	0	0.5661
No	36	31	4		6	28		31	4
AFP
+	13	13	0	0.5382	3	9	0.6639	11	1	1
−	25	21	3		4	20		23	2
Metastasis
Yes	6	6	0	1	0	6	0.3058	6	0	1
No	29	25	4		7	21		25	3
Recurrence
Yes	8	7	1	1	0	8	0.1497	8	0	0.5548
No	25	22	3		7	17		21	3
Survival
Yes	31	26	5	0.3102	7	23	0.1612	26	4	0.5558
No	10	10	0		0	10		10	0

[i] +, ΔIL>0; -, ΔIL<0. CEACAM1, carcinoembryonic antigen-related cell adhesion molecule 1.

Furthermore, the maximal information coefficient (MIC) was used to assess the dependence between clinical features and the three DS genes. We focused on three classes of splicing changes in HCC tumors, including inclusion isoform expression, IL and ΔIL. The former two only considered the degree of AS in HCC tumor samples, while ΔIL reflected the AS changes between HCC tumors and non-tumors. Then we identified their significant association with clinical features at P<0.05 (Fig. 5C). Notably, ΔIL of CEACAM1 exon 7 was significantly associated with survival time (MIC=0.38915, P<0.024; Fig. 5D). The period of overall survival time was significantly shorter in patients with increased inclusion levels of exon 7 in HCC samples (P<0.014; Fig. 5E). Moreover, ΔIL of VPS29 exon 2 was closely associated with overall survival time (MIC=0.38142, P<0.029; Fig. 5F) and cell differentiation stages (MIC=0.47849, P<0.002; Fig. 5G).

Discussion

AS plays a critical role in both physiological and pathological processes. The dysregulation of AS is highly associated with diseases such as cancer and neurodegenerative diseases (26). Recently, significant progress has been made in deciphering the role of aberrant AS events in carcinogenesis and drug resistance, and several novel molecules have been identified as possible therapeutic targets in clinical HCC research. Splicing-based prognostic biomarkers as well as therapeutic options hold great potential towards improvements in cancer therapy (27).

In the present study, we examined the relationship between AS events and human HCC and described several aberrant AS events in clinical HCC samples, such as increased inclusion of exon 2 of VPS29 and the exclusion of exon 3 of ISOC2, which could aid in HCC diagnosis and therapy. One important aberrant AS event existed in CEACAM1, exon 7 which revealed increased inclusion levels in HCC samples. CEACAM1 is a member of the immunoglobulin superfamily and a type-I transmembrane glycoprotein of the carcinoembryonic antigen family (28). In HCC, the loss of CEACAM1 indicates high metastatic potential and is associated with poor survival even after hepatectomy (29). However, research on the genome scale indicated that the different CEACAM1 isoforms had divergent roles during tumorigenesis. CEACAM1-L (CEACAM1 long cytoplasmic domain isoform) is a risk factor for HCC recurrence and its overexpression is involved in tumor cell invasion and metastasis (30).

More specifically, in this study, we demonstrated that exon 7 of CEACAM1 was the critical hub of the functional differences between the CEACAM1 isoforms. Patients with increased inclusion levels of CEACAM1 exon 7 in HCC samples had significantly shorter survival times. Exon 7 is located on the 3′terminal of the CEACAM1 transcript, although it does not encode the immunoglobulin-like domain, the inclusion or exclusion of which may influence the subcellular location of the protein, which could remodel the tumorigenicity of CEACAM1 isoforms (30–32).

Previous association analyses between candidate biomarkers and clinical features mainly utilized contingency tables and Fisher's exact test (15,33,34). Some qualitative features, such as sex, viral status and cell differentiation could not well observe the potential relevant trends with expression or splicing changes. In the present study, the MIC was first applied to the clinicopathological analysis, including both qualitative and quantitative features. This MIC method efficiently characterized the relationships for all clinical features and is expected to extend the applications in the biomarker research field.

In summary, the present study undertook a strategy of ‘small sample sequencing and large sample validation’, which means that a small number of samples for deep sequencing provided the differential spliced biomarker candidates, and a large number of samples for validation experiments confirmed the significantly differential splicing. The development of high-throughput technology has led to the generation of huge amounts of data on HCC and other cancer types. We anticipate that our method could be widely applied to thoroughly research these big datasets and discover more valuable biomarkers for HCC. The novel splicing biomarkers in this study could also provide a new understanding of HCC research.

Acknowledgements

Not applicable.

Funding

The present study was partially supported by the National Key R&D Program of China (nos. 2016YFC0902400 and 2017YFC0906603), the National Natural Science Foundation of China (nos. 81770581, 81570526 and 81123001), the Innovation project (16CXZ027), the Program of International S&T Cooperation (nos. 2014DFB30020 and 2014DFB30010), the Beijing Science and Technology Project (Z161100002616036), the Open Project Program of the State Key Laboratory of Proteomics (Academy of Military Medical Sciences, SKLP-O201509). All funding was received during this study, and there was no additional external funding received for this study.

Availability of data and materials

The datasets used during the present study are available from the corresponding author upon reasonable request.

Authors' contributions

PW and DZ designed the study, analyzed the results and created a draft of the manuscript; YW and WL collected the experimental samples and performed the qPCR experiment; AS, HW and YF prepared the experimental materials and performed the statistical analysis; XC and YJ conceived the study, performed the data interpretation, revised and approved the manuscript. All authors read and approved the manuscript and agree to be accountable for all aspects of the research in ensuring that the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Ethics approval and consent to participate

All experimental protocols were approved by the Clinical Research Ethics Committee of China-Japan Union Hospital, Jilin University.

Patient consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Glossary

Abbreviations

Abbreviations:

References