CBX8 and CD96 Are Important Prognostic Biomarkers of Colorectal Cancer.

BACKGROUND Colorectal cancer (CRC) is one of the most common malignancies worldwide, with high morbidity and mortality rates. The purpose of this study was to identify potential biomarkers in the progression of CRC. MATERIAL AND METHODS Gene and isoform expression datasets of CRC was downloaded from The Cancer Genome Atlas (TCGA). EBSeq of R was used for the normalization of gene and isoform expression, as well as the identification of differential expression genes (DEGs) and isoforms (DEIs) of CRC samples compared with normal samples. The enriched functions of DEGs and DEIs were obtained based on the Database for Annotation, Visualization and Integrated Discovery (DAVID). An independent dataset, GSE38832, was downloaded from the Gene Expression Omnibus (GEO) database for survival analysis of genes with sustained decreased/increased expression values at both gene and isoform levels with the development of CRC. RESULTS A total of 2301 genes and 4241 isoforms were found to be significantly differentially expressed in stage I-IV CRC samples. They are closely associated with muscle or cell system activity. Sixteen genes were screened out with sustained decreased/increased expression values at both gene and isoform levels with the development of CRC. Aberrant CBX8 and CD96 expressions were found to be significantly associated with CRC survival. CONCLUSIONS Through combined analysis of gene and isoform expression profiles, we identified several potential biomarkers that may play an important role in the development of CRC and could be helpful in its early diagnosis and treatment.

Background

Colorectal cancer (CRC) is one of leading causes of cancer-related death worldwide, particular in developed countries, such as the United States. Approximately 1.2 million new cases are diagnosed every year [1]. Many factors were thought to be associated with its progression and recurrence, such as stromal fibroblasts and macrophages, as well as race [2,3]. Age was also identified as a valuable predictor of prognosis and it would be better to use different strategies for patients of different ages [4,5]. Early detection was reported to effectively reduce CRC mortality and improve prognosis [6].

CRC is a heterogeneous disease whose progression is associated with multiple factors. It has been studied in many genetic and epigenetic studies and several valuable biomarkers that might contribute its initiation and recurrence have been identified. Estrogen receptor (ER) is an important transcription factor which has been verified to be implicated in cancers in many studies. Williams et al. reported that overexpression of ERβ promotes progression of CRC, and ERβ was considered as a new chemopreventive target [7]. Deregulation of gene expression in cancers are affected by many factors, including the post-transcriptional program [8]. For example, alternative splicing of Rad51C can result in overexpression of its isoform and induce the progression of CRC [9]. One of the isoforms of LD1 – LD1B – which is generated by alternative splicing, plays an important role in the maintenance of cell proliferation, and its deregulation is associated with acquisition of stem-like properties of cancer cells [10]. In addition, the post-transcriptional program can change the expression of isoforms without gene-level expression changes [8]. Therefore, it is important to explore changes at both isoform and gene levels in cancers, which has rarely been conducted.

Rapid development of gene microarray and next-generation sequencing (NGS) provides unprecedented opportunities to study cancers at the gene and isoform levels. The Affymetrix exon array contains ~5 million probes spanning ~1.4 million transcripts. It can be used to detect expression profiles at gene and exon levels based on the development of experimental and bioinformatics methods [11,12]. Some potential biomarkers in cancers have been identified based on exon microarray or NGS. For example, through the combination of exon microarray and RNA sequencing (RNA-Seq), Hoff et al. found that the fusion of VWA2-TCF7L2, DHX35-BPIFA2, and CASZ1-MASP2, as well as some novel transcript structures, may play an important role in CRC [13]. Compared with NGS, microarray has some obvious shortcomings, including the fact that it cannot detect unknown alterations of genomic structure that might contribute to progression of cancers [14]. Therefore, it may be better to use NGS in cancer research.

The Cancer Genome Atlas (TCGA, ), which was developed by the National Cancer Institute and the National Human Genome Research Institute, manipulates multiple tumor types and their corresponding genomic landscape, including gene expression profiles, as well as genomic structure variation [15]. TCGA facilitates the development of the understanding of mechanisms underlying the progression of cancers, which is limited by the small sample size. Detailed clinical information on all cancer samples is also provided, which allows the determination of associations between genes/pathways and cancer progression. In the present study, through combined analysis of isoform and gene expression profiles of CRC from TCGA, we sought to identify biological processes and genes involved in the development of CRC, which should be helpful for its early diagnosis and precise therapy.

Material and Methods

Gene and isoform expression datasets

Gene and isoform expression datasets in this study were downloaded from TCGA on 12 June 2015. We used the CRC level 3 RNA-Seq V2 datasets, which includes gene- and isoform-level expression values based on Illumina HiSeq 2000. Clinical data from CRC samples was also obtained. The normal samples were separated from CRC samples based on the TCGA barcode, and CRC samples in different stages were also separated from each other.

Differential expression analysis

Differential expression analysis and comparison of both genes and isoforms between CRC and the normal samples was conducted using empirical Bayes hierarchical models based on EBSeq [16]. Genes and isoforms with fold change >2 and FDR adjusted P-value <0.05 were screened out for the following analysis.

Functional enrichment analysis

To explore potential functions of differential expression genes (DEGs) and isoforms (DEIs), we conducted functional enrichment analysis using the Database for Annotation, Visualization, and Integrated analysis (DAVID, ) [17]. Gene Ontology (GO) terms and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways with P-value <0.05 were considered as significantly enriched.

Survival analysis

In this study, genes with sustained decreased/increased expression values with the development of CRC at both gene and isoform level were considered to be important. To explore their relations with CRC prognosis, we downloaded an independent CRC dataset from the Gene Expression Omnibus (GEO) database with the accession number GSE38832, which was deposited by Tripathi et al. [18] and consists of 122 CRC samples. Samples were divided into 2 groups based on the median expression value of a specific gene. Kaplan-Meier (KM) analysis was performed for relations of specific gene expression with CRC overall survival (OS) and the log-rank test was used for estimation of significance.

Real-time fluorescence quantitative PCR (RT-PCR)

For RT-PCR analysis, total RNA of CRC cells HCT116, which was adherent with an epithelial morphology, and normal colorectal cells NCM460, which expressed colonic epithelial cell associated antigens, was isolated with TRIzol reagent (Invitrogen, California, USA). The cDNA was synthesized with Moloney murine leukemia virus (M-MLV) reverse transcriptase (Roche, Basel, Switzerland) using random primers and was subjected to PCR amplification using rTaq polymerase (Takara Bio Inc, Shiga, Japan).

The 20-μL PCR reaction mixture contained the following: 1 μL DNA (100 ng), 0.4 μL forward primer (10 μM), 0.4 μL reverse primer (10 μM), and 10 μL SYBR Premix ExTapTM II (2×), and 8.2 μL RNase-free H2O (Roche, Basel, Switzerland).

Results

Normal and cancer samples

A total of 317 samples in TCGA with complete clinical information were obtained, which include 41 normal and 276 CRC samples. Among the CRC samples, there are 45 stage I, 110 stage II, 82 stage III, and 39 stage IV.

Differential expression genes and isoforms

DEGs and DEIs were screened out for stage I, stage II, stage III, and stage IV CRC samples compared with normal samples. Overlapping gene and isoform numbers among the 4 lists of DEGs and DEIs are shown in Figure 1. The hierarchical clustering of the overlapping DEGs/DEIs and samples is shown in Figure 2. Although the cancer samples of different stages were mixed with each other, the normal samples were mostly separated from the cancers based on the expression profiles of both DEGs and DEIs.

Figure 1

The overlapping DEGs (A) and DEIs (B) among the 4 stages of CRC samples compared with the normal ones.

Figure 2

The hierarchical clustering. (A) The hierarchical clustering of DEGs and the normal and CRC samples. (B) The hierarchical clustering of DEIs and the normal and CRC samples.

Enriched functions

A total of 49 and 29 GO terms were found to be significantly enriched in the overlapping genes and isoforms, respectively, which are mainly associated with cell and muscle system processes, as well as RNA synthesis. The top 10 most significantly enriched GO terms of the overlapping DEGs and DEIs are shown in Tables 1 and 2, respectively. DEGs and DEIs were all associated with muscle development, while RNA expression regulation-related GO terms were specific to DEGs and adhesion-related GO terms were specific to DEIs. A KEGG pathway (fatty acid metabolism, P-value=0.378) was found to be significantly enriched in overlapping DEGs, and no KEGG pathway was significantly associated with overlapping DEIs.

Table 1

The top 10 most significant GO terms of the overlapping DEGs.

Category	GO Name	P-value	Genes
BP	Regulation of transcription from RNA polymerase II promoter	6.25×10−3	MDFI, KAT2A, CBX4, NUFIP1, LEF1, PDX1, PPARGC1A, BRCA1, TARBP1, LIF, INHBA, GTF2F2, NKX2-5, TWIST1, ALX1
BP	Regionalization	9.20×10−3	MDFI, KAT2A, HOXC13, TBX20, LEF1, EN1, ALX1
BP	Regulation of striated muscle tissue development	0.0105	BCL2, LEF1, NKX2-5, TWIST1
BP	Chordate embryonic development	0.0107	KAT2A, TBX15, LEF1, EN1, HOXD1, NKX2-5, BRCA1, TWIST1, ALX1
BP	Pattern specification process	0.0110	MDFI, KAT2A, HOXC13, TBX20, LEF1, EN1, GRHL3, ALX1
BP	Regulation of muscle development	0.0111	BCL2, LEF1, NKX2-5, TWIST1
BP	Embryonic development ending in birth or egg hatching	0.0112	KAT2A, TBX15, LEF1, EN1, HOXD1, NKX2-5, BRCA1, TWIST1, ALX1
BP	Positive regulation of transcription, DNA-dependent	0.0116	LIF, INHBA, GTF2F2, TBX20, NUFIP1, LEF1, PDX1, PPARGC1A, NKX2-5, BRCA1, ALX1
MF	Neurotransmitter receptor activity	0.0119	GABRD, ANXA9, SSTR1, CHRNA6, NPFFR1
BP	Positive regulation of RNA metabolic process	0.0123	LIF, INHBA, GTF2F2, TBX20, NUFIP1, LEF1, PDX1, PPARGC1A, NKX2-5, BRCA1, ALX1

BP – biological process; MF – molecular function.

Table 2

The top 10 most significant GO terms of the overlapping DEIs.

Category	GO Name	P-value	Genes
BP	Muscle contraction	2.32×10−4	LTB4R, MYL6B, CALD1, TAZ, ASPH, HOMER1, FKBP1B
BP	Muscle system process	3.84×10−4	LTB4R, MYL6B, CALD1, TAZ, ASPH, HOMER1, FKBP1B
BP	Cell adhesion	5.57×10−4	TYRO3, PDPN, KITLG, MFGE8, CCR8, CD96, COL17A1, CORO1A, PTK2B, ITGAV, COMP, COL8A1, ADAM15
BP	Biological adhesion	5.64×10−4	TYRO3, PDPN, KITLG, MFGE8, CCR8, CD96, COL17A1, CORO1A, PTK2B, ITGAV, COMP, COL8A1, ADAM15
BP	Cell proliferation	8.00×10−4	PDPN, PTK2B, VEGFA, CD276, TNFSF14, KITLG, COL8A1, TACC3, FKBP1B, ERCC2
CC	Basement membrane	0.00171	COL17A1, VEGFA, COL8A1, ENTPD2, RELL2
BP	Cell-substrate adhesion	0.00234	CORO1A, COL17A1, PTK2B, ITGAV, ADAM15
CC	Extracellular matrix part	0.00733	COL17A1, VEGFA, COL8A1, ENTPD2, RELL2
CC	Postsynaptic density	0.0114	KLHL17, PTK2B, CALD1, HOMER1
BP	Positive regulation of cell migration	0.0142	CORO1A, PDPN, PTK2B, VEGFA

BP – biological process; CC – cellular component.

Prognostic biomarkers for CRC development and prognosis

A total of 143 DEIs and 192 DEGs (data not shown) were found to be steadily down- or up-regulated with the development of CRC (from stage I to stage IV). In addition, 16 genes were found to be steadily down- or up-regulated at both gene and isoform levels (down-regulated genes: C20orf151, CD96, ZC3H12D, MOSC2 and IER5L; up-regulated genes: NAT1, LOC388796, EHHADH, LPCAT3, DNAJC28, PHF7, KBTBD11, WDR90, MSH5, CBX8, and FXYD5). As an example, the expression values of CBX8 and CD96 at gene and isoform levels in normal and cancer samples of different stages are shown in Figure 3. The expression values of CBX8 increased with the progression of CRC at gene (Figure 3A) and isoform (Figure 3B) level. For CD96, the expression values at gene (Figure 3C) and isoform (Figure 3D) levels decreased with CRC development. High CBX8 (Figure 4A) and CD96 (Figure 4B) levels were found to be significantly associated with worse and better CRC OS, respectively.

Figure 3

The expression value of CBX8 and CD96 at gene and isoform levels in normal samples and in different stages of CRC. (A) Expression value of CBX8 at gene level. (B) Expression value of CBX8 at isoform level. (C) Expression value of CD96 at gene level. (D) Expression value of CD96 at isoform level.

Figure 4

The KM survival plot of CBX8 (A) and CD96 (B) obtained from GSE38832.

RT-PCR analysis

For CBX8 and CD96, we further confirmed their expression differences in CRC cells HCT116 and colorectal cells NCM460 through use of RT-PCR. Primers for CBX8 and CD96 were as follows: CBX8: Forward Primer: ATGGAGCTTTCAGCGGTGG, Reverse Primer: ATGCGTCCTTTCCGTATGCG; CD96: Forward Primer: CAAACACAGACAGTAGGCTTCTT Reverse Primer: GGGGATGATAGACAGCAATCAG. As a result, consistent with the results of RNA-seq analysis, the expression level of CBX8 in HCT116 cells was significantly higher than in NCM460 cells, and expression level of CD96 in HCT116 cells was significantly lower than in NCM460 cells (Figure 5).

Figure 5

RT-PCR analysis of CBX8 (A) and CD96 (B) expression levels in CRC cell line HCT116 and colorectal cell line NCM460.

Discussion

Rapid development of high-throughput gene expression profiling technologies, such as gene microarray and NGS, have accumulated a large volume of transcriptome data, most of which remains unexploited. Some databases have been developed to deposit them and can be accessed freely, such as GEO, Sequence Read Archive (SRA), and TCGA. Deep analysis of these data would revolutionize understanding of the progression of many diseases, particularly cancers. Here, through analysis of the many CRC samples from TCGA, we sought to identify some potential biomarkers that play important roles in CRC.

CRC is one of the most common malignant cancers, with ~20% distant metastasis and a low 5-year survival rate [19]. CRC initiation can attribute chromosomal instability, aberrant DNA methylation, and microsatellite instability [20]. Inflammation and microRNAs were also reported to be associated with CRC prognosis and survival [21]. Some risk factors that contribute to the progression of CRC have been identified, including diabetes, body mass index (BMI), and intestinal microenvironment [22,23]. Genetics and epigenetics were also verified to contribute to the initiation and progression of CRC. For example, down-regulation of RAB27A and RAB27B is associated with metastasis and poor prognosis [24], and 2 single-nucleotide polymorphisms (SNP) – rs2234767 and rs1800682 – of FAS can affect the progression of CRC by regulating the recruitment to chromatin of SP1/STAT1 complex [25]. Strikingly, muscle and skeletal system-related processes were found to be significantly enriched in both DEGs and DEIs in our study. Previous studies have proven skeletal and muscle metastasis sites in multiple types of cancers [26,27]. Variations of some muscle-specific molecules were also verified to affect the progression of CRC. For example, the muscle-specific miR-113b was down-regulated in CRC [28], and the expression of alpha smooth muscle actin (caspase-3 and Ki-67) can influence the effectiveness of Ocoxin® oral solution in liver metastatic CRC [29]. All of those studies should provide valuable therapeutic targets for CRC.

Alternative splicing is an important post-transcriptional mechanism that can result in multiple isoforms of a single gene. Deregulation of an isoform sometimes has no influence on the overall expression of the mRNA, but still facilitates the progression of some diseases. In addition, rapid development of NGS promotes the identification of more and more new isoforms, some of which are thought of potential biomarkers in cancers. For example, the ELF5 isoform is tissue-specific and is significantly altered in many types of cancers [30], and the alternative splicing of ER, HER2, and CD44 were also screened out in a study by Inoue [31]. In the present study, we identified 16 genes with sustained down- or up-regulated expression values at both isoform and gene levels with the development of CRC, which may be important biomarkers of CRC. Among these 16 genes, expression levels of CBX8 and CD96 were illustrated in normal and in cancer samples in different stages (Figure 3). Combined with the survival analysis from an independent dataset (Figure 4), our results suggest that CBX8 acts as an oncogene and CD96 acts as a suppressor gene in CRC.

CBX8, also known as PC3 and RC1, is a novel biomarker with DNA repair function and can promote tumorigenesis in several types of cancers [32]. CBX8 was found to be overexpressed in hepatocellular carcinoma (HCC) and was significantly associated with poor prognosis of HCC patients [33]. In CRC, the role of CBX8 remains unclear due to its function of inhibiting CRC cell proliferation but promoting CRC metastasis [34]. We identified its overexpression in CRC and its relationship with poor CRC prognosis in different datasets for the first time. To determine the role of CBX8 in CRC, further experimental studies are still needed. For CD96 (TACTILE), studies have mainly focussed on acute myeloid leukemia [35-37], but it might be a novel target in the development of CRC.

Conclusions

Through combined analysis of expression data at gene and isoform levels, we screened out 16 signatures with sustained decreased/increased expression values with CRC development. Survival analysis based on another independent dataset also indicated important roles of some of the 16 genes in CRC. Our study provides some potential targets that may be helpful in the treatment of CRC.