This article is mentioned in:

Introduction

Alternative splicing of precursor mRNAs (pre-mRNAs) is critical for regulating transcriptome diversity and protein multiplicity. During the process of alternative splicing, exons are joined together to form different transcripts, leading to the synthesis of many more proteins (1). It has been shown that alternative splicing occurs widely in eukaryotes. For example, 95% of intron-containing genes undergo alternative splicing in humans (2). There are generally five different mechanisms of alternative splicing: exon skipping, mutually exclusive exons, alternative 3′ splice site, alternative 5′ splice site and intron retention. Among them, exon skipping is the most common mode, that is, almost 35% of human alternative splicing is caused by exon skipping. In addition to those basic modes of alternative splicing, there are some other methods such as multiple promoters and multiple polyadenylation sites in eukaryotes that also occur (3–5).

Abnormally spliced mRNAs are found in many diseases, particularly in cancers. The number and types of alternative splicing differ in cancer cells compared with normal cells. For instance, cancer cells exhibit higher levels of intron retention but lower levels of exon skipping (6–9).

Colon cancer is the third leading cause of cancer-related deaths worldwide. Several studies have investigated the role of alternative splicing in colon cancer progression, and discovered that the number of alternative splicing increases during the transition from normal colon tissue to primary tumor, but decreases during metastasis to the liver (10–13). The differentially expressed alternative splicing genes in colon cancer are involved in cell-cell and cell-matrix interactions (14). Also, some alternative splicing genes have been identified because of their close association with cell growth and invasion in colon cancer, including SLC39A14, VEGF, CyclinD1, VCL, CALD1, and B3GNT6 (14,15). Since alternative splicing is crucial in colon cancer development, it is a promising target for novel anti-colon cancer therapeutics.

Preoperative chemoradiation therapy is a common and effective approach for cancer therapy, especially in colon cancer (16). Many complicated cellular responses are involved in chemoradiation therapy, which induce cancer cell death (17). Chemoradiation resistance that develops during treatment may be caused by several genetic aberrations, such as a p53 mutation and thymidylate synthase overexpression (18) However, the function of alternative splicing events in chemoradiation resistance remains unclear. In this study, we implemented a genome-wide transcriptome sequencing in HCT116 and chemoradiation resistant HCT116 (RCR-HCT116) cells to identify the alternative splicing events that affect tumor sensitivity to chemoradiotherapy.

Materials and methods

Cell culture and the construction of the RCR-HCT116 cells

HCT116 cells were cultured in McCoy's 5A (Sigma-Aldrich, St. Louis, MO, USA) containing 10% fetal bovine serum (FBS) (Gibco, Carlsbad, CA, USA) at 37°C under 5% CO2 in a CO2 incubator (Thermo Labsystems, Vantaa, Finland). HCT116 cells were exposed to 6 MV X-rays (4 Gy) at room temperature, followed by treatment with 10 µM 5-fluorouracil (5-FU) for 24 h to induce the apoptosis of tumor cells. Then the medium was replaced with fresh medium and cultured until the cells were recovered. The cells were treated with the same aforementioned method, nine times. The chemoradiation-treated cells were passaged and expanded to generate the RCR-HCT116 cells.

RNA extraction, library construction, and sequencing

Total RNA from the HCT116 and RCR-HCT116 cells were extracted using RNeasy Mini kit (Qiagen, Valencia, CA, USA) according to the manufacturer's protocol, and treated with RNase-free Dnase I (Invitrogen Life Technologies, Carlsbad, CA, USA) to remove genomic DNA.

A total amount of 2 µg RNA per sample was used for the construction of cDNA libraries. The cDNA libraries were generated according to the TruSeq® RNA Sample Prep kit v2 (Illumina, San Diego, CA, USA) following the manufacturer's protocol and then sequenced by Illumina HiSeq 4000 platform (Illumina) to generate 150-bp paired-end reads.

Read filtering and mapping

RNA-seq raw data of the HCT116 and RCR-HCT116 cells were cleaned by removing the adaptor sequence and low quality reads (mapping quality <20). Clean reads were aligned to the human reference genome sequence hg19 using TopHat (19). The following parameters were set: maximum number of mismatches permitted, 2; maximum alignments allowed, 20; maximum number of mismatches permitted in each segment alignment for reads mapped independently, 2; maximum insertion length, 3; maximum deletion length, 3; maximum mismatches in the anchor region, 0; min isoform-fraction, 0.15; minimum intron length, 50; maximum intron length, 50,000.

Differential expression analysis

Genes that were differentially expressed between the two groups were defined as differentially expressed genes (DEGs). Cufflinks (20) software was used to calculate the fragments per kilobase of transcript per million fragments mapped (FPKM) value of the different genes. The DEGs from the normal and RCR-HCT116 cells were identified by Cuffdiff (20,21) at a q≤0.01 and a fold change ≥2.

Identification and quantification of alternative splicing events

The alternative splicing events were classified into five patterns by the mixture of isoforms (MISO) (22), including alternative 5′ splice site, alternative 3′ splice site, mutually exclusive exon, intron retention and exon skipping. The MISO Bayesian inference model was used for the quantification of alternative splicing events. The change of splicing isoforms was analyzed using the MISO Bayesian inference model. The significant differentially spliced events were determined by Bayes' factor (BF) and Psi values (percent-spliced-in, Ψ) (BF≥10 and Ψ≥0.2).

Gene ontology (GO) analysis

The hypergeometric distribution test was used to identify GO categories (biological process, cellular compartment and molecular function) that were significantly enriched in a specified gene set. GO analysis was implemented with Go.db package (23).

Results

RNA-sequencing of normal and RCR-HCT116 cells

In order to generate chemoradiation-resistant colon cancer cells, we first exposed the colon cancer-derived HCT116 cell line to 6 MV X-rays (4 Gy), and then incubated the cells in 10 µM 5-fluorouracil (5-FU). As shown in Fig. 1A, after a 24-h incubation, the medium was replaced with fresh medium to remove the apoptotic and dead cells. The cells were cultured in fresh medium until they were recovered and then treated with the same aforementioned method nine times.

Figure 1 Flow of the experiments and data analysis. (A) The acquisition of chemoradiation-resistant HCT116 (RCR-HCT116) cells. (B) RNA-sequencing data analysis pipeline. 5-FU, 5-fluorouracil; RPKM, reads per kilobase million.

Sequencing was performed on the Illumina HiSeq 4000 platform to generate 150-bp paired-end reads. After the removal of the adaptor sequence and the reads of low quality, a total of 162,398,220 and 134,409,494 reads of 101 bp were generated from the HCT116 and RCR-HCT116 cells, respectively. There were 80.28 and 80.30% of the total reads from HCT116 and RCR-HT116 cells mapped to the human reference genome (Fig. 1B, Table I). Comparison of the RNA-seq data to the annotated human reference genome revealed that ~77% of the mapping reads were mapped to the CDS region. Meanwhile, the two RNA-seq libraries showed similar genomic distribution patterns from the mapping reads (Fig. 2A).

Figure 2 Overview of the RNA-seq data. (A) The genomic distribution pattern of the mapping RNA-sequencing reads of the HCT116 cells was similar to that of the chemoradiation-resistant RCR-HCT116 cells. 'CDS' refers to coding sequences; '5′ UTR' refers to 5′ untranslated region; '3′ UTR' refers to 3′ untranslated region; 'Intron' refers to regions mapped to the related location according to Refseq annotations; 'TSS_up_10 kb' refers to the 10 kb region upstream of the transcription start site (TSS); 'TES_down_10 kb' refers to the 10 kb region downstream of the transcription termination site (TES). (B) The relationship between sequencing depth and alternative splicing detection power. The x-axis indicates the sequencing depth. The y-axis indicates the number of identified splicing junctions. The pink lines indicate the known junctions; the green lines indicate the novel junctions.

Table I Statistics of the RNA-seq reads and mapped reads ratio against the human reference genome.

Table I

Statistics of the RNA-seq reads and mapped reads ratio against the human reference genome.

RNA-seq	Cells	Total reads	Total mapped reads	Mapped reads ratio (%)
WGC053648R	HCT116	162,398,220	130,371,020	80.28
WGC053649R	RCR-HCT116	134,409,494	107,929,116	80.30

DEG screening

To investigate the genes involved in chemoradiation resistance in colon cancer, we compared the RNA-seq data of the HCT116 and RCR-HCT116 cells to identify the DEGs. A total of 1,103 significant DEGs were identified (log2 FC>1, FPKM1+FPKM2 >1, q<0.01), including 818 genes that were lowly expressed and 285 genes that were highly expressed in the RCR-HCT116 cells (Table II).

Table II Top 10 differentially expressed genes.

Table II

Top 10 differentially expressed genes.

Downregulated genes	Log2 FC	Q-value	Upregulated genes	Log2 FC	Q-value
TACSTD2	−9.25	0	RFPL4A	7.43	6.27E-05
RPS4Y1	−7.49	0	SFTA1P	6.00	2.13E-03
PDE4B	−6.66	0	LARGE	5.57	0
SHF	−6.25	6.84E-08	KRT34	5.50	1.92E-07
HNF4A	−6.08	9.88E-06	GADD45G	5.46	1.05E-05
EHF	−5.91	0	KRTAP2-3	5.44	0
MKX	−5.83	7.87E-13	CREB5	5.35	3.72E-09
NTSR1	−5.79	0	HEPH	5.33	2.54E-09
KLHL35	−5.51	1.06E-07	SNAI2	5.32	2.04E-05
NEK3	−5.42	4.50E-09	TTYH1	5.17	1.77E-09

[i] Downregulated genes are genes that were expressed at a higher level in the HCT116 than in the RCR-HCT116 cells. Upregulated genes are genes that were expressed at a higher level in the RCR-HCT116 than in the HCT116 cells. RCR-HCT116, chemoradiation-resistant HCT116 cells.

GO enrichment analysis of DEGs

We then carried out GO enrichment analysis on the DEGs using Go.db package (23), which calculates the p-values using hypergeometric distribution. Genes that were expressed at a lower level in the RCR-HCT116 cells than the level in the normal HCT116 cells were enriched for GO categories related to cell cycle and cell division (e.g. cell cycle, cell cycle phase, M phase, cell cycle process, mitosis and nuclear division); whereas genes that were expressed at a higher level in RCR-HCT116 cells than in the normal HCT116 cells were enriched for GO categories related to the regulation of system processes, response to wounding, negative regulation of phosphorylation and regulation of phosphorylation (Table III).

Table III Gene ontology analysis on differentially expressed genes in HCT116 and RCR-HCT116 cells.

Table III

Gene ontology analysis on differentially expressed genes in HCT116 and RCR-HCT116 cells.

GO ID	Term	Count	P-value	Pop hits	Fold enrichment	Category
Go categories enriched in genes that were expressed higher in the HCT116 than the RCR-HCT116 cells
GO:0007049	Cell cycle	156	2.21E-62	776	4.56	BP
GO:0022403	Cell cycle phase	102	4.16E-48	414	5.59	BP
GO:0000279	M phase	91	2.89E-47	329	6.28	BP
GO:0022402	Cell cycle process	117	3.91E-47	565	4.70	BP
GO:0007067	Mitosis	70	1.74E-40	220	7.22	BP
GO:0000280	Nuclear division	70	1.74E-40	220	7.22	BP
GO:0048285	Organelle fission	71	3.07E-40	229	7.04	BP
GO:0000087	M phase of mitotic cell cycle	70	6.53E-40	224	7.09	BP
GO:0000278	Mitotic cell cycle	84	1.73E-36	370	5.15	BP
GO:0051301	Cell division	73	3.25E-34	295	5.62	BP
GO:0006259	DNA metabolic process	93	9.49E-33	506	4.17	BP
GO:0031981	Nuclear lumen	158	3.60E-32	1450	2.64	CC
GO:0005694	Chromosome	81	3.70E-29	460	4.26	CC
GO:0070013	Intracellular organelle lumen	172	5.91E-29	1779	2.34	CC
GO:0043233	Organelle lumen	173	2.95E-28	1820	2.30	CC
GO:0031974	Membrane-enclosed lumen	175	3.27E-28	1856	2.28	CC
GO:0043228	Non-membrane-bounded organelle	214	2.08E-27	2596	2.00	CC
GO:0043232	Intracellular non-membrane-bounded organelle	214	2.08E-27	2596	2.00	CC
GO:0000793	Condensed chromosome	43	7.82E-27	129	8.07	CC
GO:0044427	Chromosomal part	69	3.73E-25	386	4.33	CC
Go categories enriched in genes that were expressed higher in the RCR-HCT116 than the HCT116 cells
GO:0044057	Regulation of system process	18	3.18E-06	309	3.94	BP
GO:0009611	Response to wounding	22	3.58E-05	530	2.81	BP
GO:0006937	Regulation of muscle contraction	8	8.80E-05	72	7.52	BP
GO:0043005	Neuron projection	15	4.55E-04	342	3.01	CC
GO:0042326	Negative regulation of phosphorylation	6	4.96E-04	45	9.02	BP
GO:0042325	Regulation of phosphorylation	18	5.30E-04	466	2.61	BP
GO:0010563	Negative regulation of phosphorus metabolic process	6	6.71E-04	48	8.46	BP
GO:0045936	Negative regulation of phosphate metabolic process	6	6.71E-04	48	8.46	BP
GO:0051174	Regulation of phosphorus metabolic process	18	8.28E-04	485	2.51	BP
GO:0019220	Regulation of phosphate metabolic process	18	8.28E-04	485	2.51	BP
GO:0008285	Negative regulation of cell proliferation	15	9.19E-04	361	2.81	BP
GO:0042127	Regulation of cell proliferation	24	1.28E-03	787	2.06	BP
GO:0005201	Extracellular matrix structural constituent	7	1.66E-03	86	5.50	MF
GO:0042981	Regulation of apoptosis	24	1.69E-03	804	2.02	BP
GO:0043067	Regulation of programmed cell death	24	1.92E-03	812	2.00	BP
GO:0040013	Negative regulation of locomotion	6	2.00E-03	61	6.65	BP
GO:0010941	Regulation of cell death	24	2.02E-03	815	1.99	BP
GO:0040012	Regulation of locomotion	10	2.15E-03	192	3.52	BP
GO:0044459	Plasma membrane part	48	3.05E-03	2203	1.50	CC
GO:0006954	Inflammatory response	13	3.14E-03	325	2.71	BP

[i] Fold enrichment, the ratio between the predicted and observed number of genes within the given GO category. Count, the number of genes observed in the given GO category. Category: BP, 'biological process'; MF, 'molecular function'; CC, 'cellular component'; GO, Gene ontology; RCR-HCT116, chemoradiation-resistant HCT116 cells.

Identification and annotation of alternative splicing events

To determine the relationship between sequencing depth and the detection power of alternative splicing, the sequencing libraries (HCT116 and RCR-HCT116 cells) were randomly selected to create sub-libraries (i.e. 5–95% of the whole library) to determine the known and novel junctions. As shown in Fig. 2B, the sequencing depth was correlated with the detection of unknown junctions, however, when the sequencing depth was more than 20% of the whole library, the increase of sequencing depth did not significantly increase the number of known junctions. This implies that our sequencing data was capable of supporting the identification of unknown junctions.

We further examined the splicing patterns in the normal and chemoradiation-resistant colon cancer cells using the MISO (22) package. MISO quantifies the level of inclusion of a given differentially expressed exon as the 'percent spliced in' (Psi or Ψ), which reflects the fraction of a gene's mRNA that includes the exon, intron or alternative splice site. Ψ values vary between 0 (the exon, intron or alternative splice site is excluded from every transcript) and 1 (the exon, intron or alternative splice site is included in every transcript). MISO also calculates a Bayes factor for each differential splicing event, which is a measure of the odds that there is differential inclusion of a particular exon in different samples. The five main alternative splicing patterns, 3′ splice site (A3SS), alternative 5′ splice site (A5SS), mutually exclusive exon (MXE), intron retention (IR) and exon skipping (ES), were analyzed in the RNA-seq data of the HCT116 and the RCR-HCT116 cells (Fig. 3A). Fig. 3B shows the Sashimi plots of five examples with different patterns of alternative splicing events, with the number of reads that span each part of the splice junction shown on the plots for the two samples analyzed. Furthermore we calculated the number of alternative splicing events for both types of cells and as shown in Fig. 3C, there was no significant difference in the number of detected A3SS, A5SS, MXE and IR in the normal and RCR-HCT116 cells, indicating that these types of alternative splicing may not function in chemoradiation-resistant colon cancer cells. Nevertheless, the number of ES was significantly increased in the chemoradiation-resistant colon cancer cells (Fig. 3C), suggesting that chemoradiation may functions via ES.

Figure 3 Statistical analysis of the different alternative splicing events. (A) The schematic plots show the five types of alternative splicing patterns. (B) The Sashimi plots of the RNA-seq data for the five types of alternative splicing events that are expressed differentially in the HCT116 and RCR-HCT116 cells. The main panel shows the counts of RNA-seq reads that span the junctions in each region. The coordinates for each splicing event are shown at the top, and the schematic of this splicing event is shown at the bottom. Dark pink indicates the HCT116 cells, orange indicates the RCR-HCT116 cells. The estimated MISO values are shown in the right panels, the dashed lines indicate the 95% confidence intervals of the MISO values. The Sashimi plot was produced using the MISO package. (C) The number of different splicing events in the normal and chemoradiation-resistant colon cancer cells. A3SS, alternative 3′ splice site; A5SS, alternative 5′ splice site; MXE, mutually exclusive exon; IR, intron retention; ES, exon skipping; RCR-HCT116, chemoradiation-resistant HCT116 cells; MISO, mixture of isoforms; RPKM, reads per kilobase million.

Go enrichment analysis of genes with differentially alternative splicing levels

To gain further insight into the role of these alternative splicing level altered genes, we performed GO analysis on the genes that had different alternative splicing levels in the normal and RCR-HCT116 cells. Our dataset identified that 323 alternative splicing events in 293 genes were significantly different between the normal and RCR-HCT116 cells (data not shown). Some of the alternative splicing-containing genes were previously reported in colon cancers, such as CD44. Tumors carrying the CD44 v6 epitope (exon v6) acquire selective advantages during tumor progression. We further identified three new MXEs between exon 5 and exon 11 of CD44 (Fig. 4), which may be associated with chemoradiation resistance of colon cancer. Go analysis showed that these genes were clustered functionally into several groups related with DNA replication, such as deoxyribonucleotide metabolic/catabolic processes (MPG, NT5M, RRM2B, OGG1, NT5C), response to DNA damage-stimulus (POLL, RECQL4, MPG, C17ORF70, PCBP4, POLG, ZMAT3, MUS81, GTF2H4, ZSWIM7, RRM2B, OGG1, ERCC3) and helicase activity (RECQL4, MOV10, DDX11, SKIV2L, EIF4A1, SNORA67, HLTF, ERCC3, SMARCA4) (Fig. 5).

Figure 4 Differential splicing of CD44 in the HCT116 and RCR-HCT116 cells. The main panel shows the counts of RNA-seq reads that span the junctions in each region. The coordinates for each splicing event are shown at the top, and the schematic of this splicing event is shown at the bottom. Dark pink indicates the HCT116 cells and the orange indicates the RCR-HCT116 cells. The estimated MISO values are shown in the right panels, the dashed lines indicate the 95% confidence intervals of the MISO values. The Sashimi plot was produced using the MISO package. RCR-HCT116, chemoradiation-resistant HCT116 cells; MISO, mixture of isoforms RPKM, reads per kilobase million.

Figure 5 Top Gene ontology (GO) categories enriched in the genes that were found to have different alternative splicing levels in the normal and chemoradiation-resistant HCT116 cells. The black bars indicate the −log (base 10) of the p-values for the hypergeometric distribution test. The gray bars indicate the ratio between the predicted and observed number of genes within the given GO category (fold enrichment).

Discussion

Alternative splicing, a key molecular event that allows for protein diversity, is an important post-transcriptional regulatory mechanism to control cell processes. Aberrant splicing is related with various diseases, including colon cancer (24–26). The importance of alternative splicing in colon cancer progression has been emphasized in many studies (24). Until now, many alternative splicing genes have been identified in colon cancer because of their close association with cell growth and invasion, including SLC39A14, VEGF, CyclinD1, VCL, CALD1, B3GNT6, ACTN1, TPM1, FN1, COL6A3, SLC3A2 (13–15). Evidence suggests the role of CD44 alternative splicing in the progression of colon cancer. The expression of CD44 splice variants (exon v6) is increased during colon cancer progression and the expression level of CD44 v6 is associated with tumor-related mortality (27–29). In this study, we identified three novel MXEs between the exon 5 and exon 11 of the CD44 gene. The levels of these MXEs were different between the normal and RCR-HCT116 cells, indicating that these novel alternative-splicing events occurring in CD44 may be related with the chemoradiation resistance of colon cancer. Thus, different alternative splicing events, that even exist on the same gene, may have diverse functions in tumor development and therapy.

Preoperative chemoradiation therapy is increasingly used in colon cancer therapy (30,31). Some patients exhibit a marked pathologic response with standard chemoradiation treatment, however others remain non-responsive. Thus, the identification of markers that can predict sensitivity to chemoradiation is exceedingly useful to avoid unnecessary preoperative treatment. Previous studies have identified several markers that predict the sensitivity or resistance to chemoradiation therapy. A clinical study showed that thymidylate synthase (TS) is overexpressed in chemoradiation-resistant rectal cancer patients, which indicates that the level of TS in tumors is the best predictor of sensitivity to chemoradiation (31). In our study, we identifed 323 alternative splicing events in 293 genes that were significantly different between the normal and chemoradiation-resistant HCT116 cells. Notably, there were no significant differences in the expression of most of these alternative splicing affected genes (26 out of 293, data not shown). It is deducible that apart from the expression level of some crucial genes, alternative-splicing events of these genes may also affect tumor sensitivity to chemoradiotherapy.

In this study, we defined a set of 293 genes showing different alternative splicing events in a normal and chemoradiation-resistant colon cancer cell line. This group of genes were enriched in molecular functions and biological processes relevant to DNA replication, such as deoxyribonucleotide metabolic/catabolic processes and helicase activity. We identified for the first time, to the best of our knowledge, the alternative splicing events that are associated with the chemoradiation resistance of colon cancer. Thus, from a clinical point of view, our study is expected to provide insight into potential novel therapeutic targets, such as alternative splicing, to improve treatment response.

Acknowledgments

This study was supported by the Joint Funds of the Department of Science and Technology of Yunnan Province and Kunming Medical University (no. 2013FB167) and the Key Project of the Department of Education of Yunnan Province (no. 2014Z061).

References