Transcriptome analysis of clock disrupted cancer cells reveals differential alternative splicing of cancer hallmarks genes

Abstract

Emerging evidence points towards a regulatory role of the circadian clock in alternative splicing (AS). Whether alterations in core-clock components may contribute to differential AS events is largely unknown. To address this, we carried out a computational analysis on recently generated time-series RNA-seq datasets from three core-clock knockout (KO) genes (ARNTL, NR1D1, PER2) and WT of a colorectal cancer (CRC) cell line, and time-series RNA-seq datasets for additional CRC and Hodgkin's lymphoma (HL) cells, murine WT, Arntl KO, and Nr1d1/2 KO, and murine SCN WT tissue. The deletion of individual core-clock genes resulted in the loss of circadian expression in crucial spliceosome components such as SF3A1 (in ARNTL^KO), SNW1 (in NR1D1^KO), and HNRNPC (in PER2^KO), which led to a differential pattern of KO-specific AS events. All HCT116^KO cells showed a rhythmicity loss of a crucial spliceosome gene U2AF1, which was also not rhythmic in higher progression stage CRC and HL cancer cells. AS analysis revealed an increase in alternative first exon events specific to PER2 and NR1D1 KO in HCT116 cells, and a KO-specific change in expression and rhythmicity pattern of AS transcripts related to cancer hallmarks genes including FGFR2 in HCT116_ARNTL^KO, CD44 in HCT116_NR1D1^KO, and MET in HCT116_PER2^KO. KO-specific changes in rhythmic properties of known spliced variants of these genes (e.g. FGFR2 IIIb/FGFR2 IIIc) correlated with epithelial-mesenchymal-transition signalling. Altogether, our bioinformatic analysis highlights a role for the circadian clock in the regulation of AS, and reveals a potential impact of clock disruption in aberrant splicing in cancer hallmark genes.

Fig. 1. RNA-seq data analysis showed differences in the circadian phenotype between gene- and transcript expression values.

a Schematic representation of the pipeline for RNA-seq data analysis. The input data consisted of 222 samples from eight mammalian time-course RNA-seq data. Murine datasets from Per1/2 KO (GSE171975) and Cry1/2 KO (GSE135898) were only analysed at the gene level. b–d Classification of 24 h rhythmic features in HCT116 datasets at (b) gene level and (c) transcript level showed a decline in the total number of circadian expressed genes/transcript-sets and changes in their biotypes due to core-clock disruption. d Scatter plot represents each of the HCT116 WT and KO cells according to the total number of 24 h rhythmic genes (x-axis) and the total number of 24 h rhythmic transcripts (y-axis).

Fig. 2. Circadian expression of splicing factors and spliceosome machinery related transcripts in mammalian RNA-seq datasets.

Median normalized phase-sorted heatmaps (left) and acrophase bins (right) of SFs circadian transcripts (tr.) in (a) HCT116_WT and HCT116^KO cells (upper panel), SW480 and SW620 cell lines (lower left panel), HDMYZ and L1236 cell lines (lower right panel), b murine WT and Arntl KO from distal bronchiolar epithelium tissue (upper left panel) and from primary tracheal epithelial cells (upper right panel), murine WT and Nr1d1/2 KO from epithelial cells (lower left panel), and murine SCN WT tissue (lower right panel).

Fig. 3. Clock disruption affected rhythmicity of splicing factors in human and murine datasets.

a–c Circular plots depict the distribution of peak phases of overlapping differentially rhythmic SF transcripts between (a) HCT116_WT and HCT116_ARNTL^KO, (b) HCT116_WT and HCT116_NR1D1^KO, and (c) HCT116_WT and HCT116_PER2^KO. d–h Categorial heatmaps represent the loss of rhythmic SFs in (d) HCT116^KO vs. WT, (e) SW620 vs. SW480, (f) L1236 vs. HDMYZ, (g) mouse Arntl KO vs. WT in distal bronchiolar epithelium (left) and primary epithelial cells (right), and (h) Nr1d1/2 KO vs. WT. Green colour indicates SFs with at least one ~24 h transcript whereas grey colour indicates SFs with no rhythmic transcript. The numbers above the categorial heatmaps indicate the total number of SFs (genes) with rhythmic transcript(s).

Fig. 4. Local alternative splicing analysis revealed changes in the pattern of loss and gain of all seven splice modes in human and murine datasets.

Bar plots depict the total number of genes (x-axis) with loss and gain of splicing (a–c) in the HCT116_ARNTL^KO, HCT116_NR1D1^KO, and HCT116_PER2^KO vs. HCT116_WT cells, (d) in the SW620 vs. SW480 cells, (e) in the L1236 vs. HDMYZ cells, (f) in Arntl KO from distal bronchiolar epithelium (left) and from primary tracheal epithelial cells (right) vs. their WT, and (g) in Nr1d1/2 KO vs. WT. In each barplot, the number of genes containing transcripts with circadian expression and protein-coding biotype (dark blue), protein-coding biotype (blue), and different biotypes (light blue) is indicated.

Fig. 5. Phase-shifted spliced isoforms in HCT116 datasets are associated with hallmarks of cancer.

Differential rhythmicity analysis between rhythmic transcripts within the same cell line was carried out using DODR. a The scatter plot depicts the distribution of differentially rhythmic transcript pairs according to their phase difference and amplitude ratio in HCT116. b Decline in differentially rhythmic isoform pairs observed in all HCT116^KO. c The chord diagram represents the peak phases of transcript pairs across HCT116 under different phase shift cut-offs. d The chord diagram depicts the biotypes of phase shifted transcript pairs. e The circular plot shows the association of genes containing phase shifted transcript pairs with hallmarks of cancer.

Fig. 6. Clock disruption in HCT116 resulted in alterations in AS in genes involved in the hallmarks of cancer.

Commonly spliced and uniquely spliced candidates in HCT116^KO were mapped to the cancer hallmarks gene list. a Circular plots depict the association between spliced candidates and different cancer hallmarks (b) HNRNPM, a spliceosome machinery component lost its circadian expression in HCT116^KO cells. HCT116^KO cells also showed discrepancies in the mean expression level of HNRNPM transcripts. Expression of uniquely alternatively spliced candidate transcripts in (c) HCT116_PER2^KO, (d) HCT116_ARNTL^KO, and (e) HCT116_NR1D1^KO were plotted. Genomic region plots of MET, FGFR2 and LRRC8A transcripts represent differences in their exon composition (marked in red) compared to canonical forms. Circadian rhythmic transcripts were plotted using harmonic regression fit and arrhythmic transcripts were plotted using Loess fit in R.